From 2d9cb5cc17e66c41eb22f757dbf900533a16656e Mon Sep 17 00:00:00 2001
From: Max Herzog <max.herzog@mailbox.tu-dresden.de>
Date: Sat, 6 Jan 2024 10:30:34 +0000
Subject: [PATCH] Upload New File
---
src/TIE.jl | 374 +++++++++++++++++++++++++++++++++++++++++++++++++++++
1 file changed, 374 insertions(+)
create mode 100644 src/TIE.jl
diff --git a/src/TIE.jl b/src/TIE.jl
new file mode 100644
index 0000000..dd9a7fb
--- /dev/null
+++ b/src/TIE.jl
@@ -0,0 +1,374 @@
+module TIE
+
+using LinearAlgebra
+using SparseArrays
+using Statistics
+using FFTW
+using NLopt
+using JSON
+using NPZ
+
+
+# constants
+global const mu0 = 4*pi*1e-7
+global const qe = 1.6e-19
+global const me = 9.1095e-31
+global const hquer = 1.05459e-34
+global const c = 2.997925e8
+
+"""
+ err(p::Vector, grad::Vector)
+
+Specifies the cost function to be minimized to find the best fitting boundary condition
+TODO: Ãœbergabe der (typdefifinierten) Parameter
+"""
+function err(p::Vector, grad::Vector, Dreg, rhs0, kx2, ky2, k0, dz, Io, Iu, If, N, N2, bandlimit)
+
+ println("iteration: ", iOpt)
+ println("boundary params: ",p)
+
+ boundvalF = zeros(ComplexF64,2*(N.x+N.y+2))
+ # band limited boundary function
+ for ii=1:bandlimit
+ boundvalF[N.x+N.y+2-bandlimit+ii] = p[end+1-2*ii] - im*p[end+2-2*ii]
+ boundvalF[N.x+N.y+3+ii] = p[2*ii-1] + im*p[2*ii]
+ end
+ boundval = real.(ifft(ifftshift(boundvalF)))
+ rhs0[1:end-1,1] = boundval[1:N.x+1]
+ rhs0[end,1:end-1] = boundval[N.x+2:N.x+N.y+2]
+ rhs0[1:end-1,end] = reverse(boundval[N.x+N.y+3:2*N.x+N.y+3])
+ rhs0[1,1:end-1] = reverse(boundval[2*N.x+N.y+4:2*N.x+2*N.y+4])
+ rhs = vec(rhs0)
+
+ # solution of TIE
+ ret = Dreg \ rhs
+ phase = reshape(ret,(N2.x,N2.y))
+ ψ = sqrt.(If).*exp.(im*phase[2:end-1,2:end-1])
+
+ # duplication of domain to avoid propagation artifacts at boundary
+ ψ2 = [ψ reverse(ψ,dims=2); reverse(ψ,dims=1) reverse(ψ,dims=(1,2))]
+
+ # over- and underfocused intensities from reconstructed wave
+ IoRec2 = abs2.(ifft(fft(ψ2).*ifftshift(exp.(im*dz/(2*k0)*(kx2.^2 .+ ky2.^2)))))
+ IoRec = IoRec2[1:N.x,1:N.y]
+ IuRec2 = abs2.(ifft(fft(ψ2).*ifftshift(exp.(-im*dz/(2*k0)*(kx2.^2 .+ ky2.^2)))))
+ IuRec = IuRec2[1:N.x,1:N.y]
+
+ # deviation between reconstructed and experimental intensities
+ dev = norm(IoRec.-Io)/norm(Io)+norm(IuRec.-Iu)/norm(Iu)
+ println("deviation: ",dev)
+
+ global iOpt = iOpt+1
+
+ return dev
+end
+
+"""
+ calculate_TIE(path::String, filename::String)
+
+This is the main function of this script. It reads the declared .json-file and performs the TIE calculation
+in accordance with the settings selected in the .json-file.
+"""
+function calculate_TIE(path::String, filename::String)
+ cd(path)
+ j = open(filename) do file
+ j = JSON.parse(file)
+ end
+ # parameter
+ Ua = j["Uacc"] # acceleration voltage
+ E = me*c^2+qe*Ua # electron energy
+ k0 = sqrt(E^2-me^2*c^4)/(hquer*c) # electron wave number
+ println("acceleration voltage [kV]: ",Ua/1e3)
+
+ reg = (j["reg"])^2/k0 # reciprocal wave vector regularization parameter (2e7 for MnPtSn)
+ if reg>0
+ println("regularization length [μm]: ",1/sqrt(reg*k0)*1e6)
+ else
+ println("no regularization")
+ end
+
+ dz = j["dz"]
+ println("defocus [mm]: ",dz*1e3)
+
+ p = (x = j["px"], y = j["py"])
+
+ save = j["save"] # save result
+ poi = j["poi"] # also reconstruct Poisson solution
+
+ filename_o = j["Io"]
+ try
+ global filename_u = j["Iu"]
+ catch
+ global filename_u = ""
+ end
+ try
+ global filename_f = j["If"]
+ catch
+ global filename_f = ""
+ end
+
+ # load image data (at least 2 defocused images)
+ if ~isempty(filename_o)
+ Io = npzread(filename_o)
+ Io = Io./mean(Io)
+ N = (x = size(Io,1), y = size(Io,2))
+ if ~isempty(filename_u)
+ Iu = npzread(filename_u)
+ Iu = Iu./mean(Iu)
+ if ~isempty(filename_f)
+ If = npzread(filename_f)
+ If = If./mean(If)
+ else
+ If = (Io.+Iu)./2
+ end
+ else
+ if ~isempty(filename_f)
+ If = npzread(filename_f)
+ If = If./mean(If)
+ Iu = 2*If.-Io
+ else
+ If = mean(Io)*ones(N.x,N.y)
+ Iu = 2*If.-Io
+ end
+ end
+ end
+
+ println("image dimensions [px]: ",N.x,"x",N.y)
+
+ a = (x = N.x*p.x, y = N.y*p.y)
+ println("image dimensions [μm]: ",a.x*1e6,"x",a.y*1e6)
+
+ thick = j["thick"] #sample thickness
+
+ # grids
+ xq = 1:N.x
+ yq = 1:N.y
+ x = a.x/N.x .* ((0:N.x-1) .- N.x÷2)
+ x = reshape(x, (N.x, 1))
+ y = a.y/N.y .* ((0:N.y-1) .- N.y÷2)
+ y = reshape(y, (1, N.y))
+ kx = 2*pi*N.x/a.x.^2*x
+ ky = 2*pi*N.y/a.y.^2*y
+
+ x2 = a.x/N.x .* ((0:2*N.x-1) .- N.x) # double size for propagation / Poisson
+ x2 = reshape(x2, (2*N.x, 1))
+ y2 = a.y/N.y .* ((0:2*N.y-1) .- N.y)
+ y2 = reshape(y2, (1, 2*N.y))
+ kx2 = pi*N.x/a.x.^2*x2
+ ky2 = pi*N.y/a.y.^2*y2
+
+ ## Possion solver with pseudo Dirichlet BCs
+ if poi
+
+ av = mean(Io.+Iu)/2
+ dIdz2 = [Io .- Iu -reverse(Io .- Iu,dims=2); -reverse(Io .- Iu,dims=1) reverse(Io .- Iu,dims=(1,2))]
+ If2 = [If reverse(If,dims=2); reverse(If,dims=1) reverse(If,dims=(1,2))]
+
+ # neglect I variation solution
+ temp = ifftshift((kx2.^2 .+ky2.^2) ./((kx2.^2 .+ky2.^2).^2 .+ k0^2*reg^2) .*fftshift(fft(dIdz2./If2)))
+ if reg == 0.
+ temp[1,1] = 0.
+ end
+ temp = -k0/dz*real.(ifft(temp))
+ phasePoi1 = temp[1:N.x,1:N.y]
+
+ # const. I solution
+ temp = ifftshift((kx2.^2 .+ky2.^2) ./((kx2.^2 .+ky2.^2).^2 .+ k0^2*reg^2) .*fftshift(fft(dIdz2)))
+ if reg == 0.
+ temp[1,1] = 0.
+ end
+ temp = -k0/dz/av*real.(ifft(temp))
+ phasePoi2 = temp[1:N.x,1:N.y]
+
+ # irrotational current solution
+ temp = av*(ifftshift(kx2.^2) .*fft(temp./If2) .+ ifftshift(ky2.^2) .*fft(temp./If2))./ifftshift(kx2.^2 .+ ky2.^2)
+ temp[1,1] = 0.
+ temp = real.(ifft(temp))
+ phasePoi3 = temp[1:N.x,1:N.y]
+
+ BxPoi1 = zeros(size(phasePoi1))
+ BxPoi1[:, 2:end] = -hquer*N.x/(a.x * thick * qe) * diff(phasePoi1, dims=2)
+ BxPoi1[:, 1] = BxPoi1[:, 2]
+ ByPoi1 = zeros(size(phasePoi1))
+ ByPoi1[2:end, :] = hquer * N.y / (a.y * thick * qe) * diff(phasePoi1, dims=1)
+ ByPoi1[1, :] = ByPoi1[2, :]
+
+ BxPoi2 = zeros(size(phasePoi2))
+ BxPoi2[:, 2:end] = -hquer*N.x/(a.x * thick * qe) * diff(phasePoi2, dims=2)
+ BxPoi2[:, 1] = BxPoi2[:, 2]
+ ByPoi2 = zeros(size(phasePoi2))
+ ByPoi2[2:end, :] = hquer * N.y / (a.y * thick * qe) * diff(phasePoi2, dims=1)
+ ByPoi2[1, :] = ByPoi2[2, :]
+
+ BxPoi3 = zeros(size(phasePoi3))
+ BxPoi3[:, 2:end] = -hquer*N.x/(a.x * thick * qe) * diff(phasePoi3, dims=2)
+ BxPoi3[:, 1] = BxPoi3[:, 2]
+ ByPoi3 = zeros(size(phasePoi3))
+ ByPoi3[2:end, :] = hquer * N.y / (a.y * thick * qe) * diff(phasePoi3, dims=1)
+ ByPoi3[1, :] = ByPoi3[2, :]
+
+ if save
+ npzwrite("phase_gradIconst.npy", phasePoi1)
+ npzwrite("phase_Iconst.npy", phasePoi2)
+ npzwrite("phase_irr.npy", phasePoi3)
+ npzwrite("B_gradIconst.npy", BxPoi1, ByPoi1)
+ npzwrite("B_Iconst.npy", BxPoi2, ByPoi2)
+ npzwrite("B_irr.npy", BxPoi3, ByPoi3)
+ end
+
+ end
+
+ ## discretized TIE with Dirichlet BCs ##
+ # sampling points for augmented image containing boundary values
+ N2 = (x = N.x+2, y = N.y+2)
+
+ # Intensity
+ temp = zeros(Float64,N2.x,N2.y)
+ temp[2:end-1,2:end-1] .= If
+ I = spdiagm(0 => vec(temp)) # sparse matrix trafo
+
+ # x-derivatives
+ temp0 = ones(Float64, N2.x*N2.y)
+ temp1 = -ones(Float64, N2.x*N2.y-1)
+ Dx1 = N.x/a.x*spdiagm(-1 => temp1, 0 => temp0)
+ Dx2 = N.x/a.x*spdiagm(0 => -temp0, 1 => -temp1)
+
+ # y-derivatives
+ temp1 = -ones(Float64, N2.x*N2.y-N2.x)
+ Dy1 = N.y/a.y*spdiagm(-N2.x => temp1, 0 => temp0)
+ Dy2 = N.y/a.y*spdiagm(0 => -temp0, N2.x => -temp1)
+
+ # differential operator
+ D = 1/2 .*(Dx2*I*Dx1 .+ Dy2*I*Dy1 .+ Dx1*I*Dx2 .+ Dy1*I*Dy2)
+
+ # regularization
+ if reg>0
+ temp = ones(Float64, N2.x*N2.y)
+ Dreg = (transpose(D)*D .+ reg^2*k0^2*spdiagm(0 => temp))
+ else
+ Dreg = D
+ end
+
+ # Dirichlet BCs
+ for i=1:N2.x
+ Dreg[i,:] .= 0.
+ Dreg[i,i] = 1.
+ end
+ for i=N2.x*N2.y-N2.x+1:N2.x*N2.y
+ Dreg[i,:] .= 0.
+ Dreg[i,i] = 1.
+ end
+ for i=1:N2.x:N2.x*N2.y
+ Dreg[i,:] .= 0.
+ Dreg[i,i] = 1.
+ Dreg[i+N2.x-1,:] .= 0.
+ Dreg[i+N2.x-1,i+N2.x-1] = 1.
+ end
+
+ dropzeros!(Dreg)
+
+ # right hand side (inhomogeneous term) of TIE
+ rhs0 = zeros(Float64,N2.x,N2.y)
+ rhs0[2:end-1,2:end-1] .= k0*(Io.-Iu)./(2*dz)
+ if reg>0
+ rhs0 = reshape(transpose(D)*vec(rhs0),(N2.x,N2.y))
+ end
+
+ # 1st run to determine bounds of boundary fourier coefficients
+ bandlimit = j["bandlim"] # number of non-zero fourier coefficient of boundary function
+ println("band limit: ",bandlimit)
+
+ boundvalF = zeros(ComplexF64,2*(N.x+N.y+2)) # boundary function set to zero
+ boundval = real.(ifft(ifftshift(boundvalF)))
+ rhs0[1:end-1,1] = boundval[1:N.x+1]
+ rhs0[end,1:end-1] = boundval[N.x+2:N.x+N.y+2]
+ rhs0[1:end-1,end] = reverse(boundval[N.x+N.y+3:2*N.x+N.y+3])
+ rhs0[1,1:end-1] = reverse(boundval[2*N.x+N.y+4:2*N.x+2*N.y+4])
+ rhs = vec(rhs0)
+
+ # solution of TIE
+ ret = Dreg \ rhs
+ phase = reshape(ret,(N2.x,N2.y))
+
+ # approximation of upper bounds for fourier coefficients of boundary function from parseval's theorem
+ bounds = 2*sqrt((maximum(phase)-minimum(phase))^2*(2*(N.x+N.y+2))^2/2)
+ println("Fourier component bounds: ",bounds)
+
+ ## optimization of boundary conditions ##
+ # global optimizers: :GN_DIRECT_L, GD_STOGO, GN_CRS2_LM
+ # local optimizers: :LN_COBYLA, :LN_SBPLX, :LN_NELDERMEAD, :LN_BOBYQA
+
+ opt = Opt(:GN_DIRECT_L, 2 * bandlimit)
+ opt.lower_bounds = -bounds .*ones(Float64, 2*bandlimit)
+ opt.upper_bounds = bounds .*ones(Float64, 2*bandlimit)
+
+ #opt.ftol_abs = 1e-15
+ opt.xtol_rel = 1e-15
+ opt.maxeval = j["maxiter"]
+ println("maximal number of iterations: ",opt.maxeval)
+
+ opt.min_objective = ((x,grad) -> err(x, grad, Dreg, rhs0, kx2, ky2, k0, dz, Io, Iu, If, N, N2, bandlimit))
+ #p0 = zeros(Float64, 2*bandlimit) # starting values for boundary function fourier coefficients
+ p0 = rand(Float64, 2*bandlimit)
+ global iOpt = 1
+
+ (minf,minx,retu) = optimize(opt, p0)
+ # global p0org=minx
+ numevals = opt.numevals # the number of function evaluations
+ println("got $minf at $minx after $numevals iterations (returned $retu)")
+ #ftbval = minx
+ if (maximum(abs.(minx))-bounds)/bounds>0.95
+ println("warning: bounds for parameter search too small")
+ end
+
+ boundvalF = zeros(ComplexF64,2*(N.x+N.y+2))
+ # band limited boundary function
+ for ii=1:bandlimit
+ boundvalF[N.x+N.y+2-bandlimit+ii] = minx[end+1-2*ii] - im*minx[end+2-2*ii]
+ boundvalF[N.x+N.y+3+ii] = minx[2*ii-1] + im*minx[2*ii]
+ end
+ boundval = real.(ifft(ifftshift(boundvalF)))
+
+ rhs0[1:end-1,1] = boundval[1:N.x+1]
+ rhs0[end,1:end-1] = boundval[N.x+2:N.x+N.y+2]
+ rhs0[1:end-1,end] = reverse(boundval[N.x+N.y+3:2*N.x+N.y+3])
+ rhs0[1,1:end-1] = reverse(boundval[2*N.x+N.y+4:2*N.x+2*N.y+4])
+
+ rhs = vec(rhs0)
+
+ # solution of TIE
+ ret = Dreg \ rhs
+ phase = reshape(ret,(N2.x,N2.y))
+
+ # final deviation between reconstructed and experimental intensities
+ ψ = sqrt.(If).*exp.(im*phase[2:end-1,2:end-1])
+ ψ2 = [ψ reverse(ψ,dims=2); reverse(ψ,dims=1) reverse(ψ,dims=(1,2))]
+ Iorec2 = abs2.(ifft(fft(ψ2).*ifftshift(exp.(im*dz/(2*k0)*(kx2.^2 .+ ky2.^2)))))
+ Iorec = Iorec2[1:N.x,1:N.y]
+ Iurec2 = abs2.(ifft(fft(ψ2).*ifftshift(exp.(-im*dz/(2*k0)*(kx2.^2 .+ ky2.^2)))))
+ Iurec = Iurec2[1:N.x,1:N.y]
+
+ dev = norm(Iorec[2:N.x-1,2:N.y-1].-Io[2:N.x-1,2:N.y-1])/norm(Io[2:N.x-1,2:N.y-1])+norm(Iurec[2:N.x-1,2:N.y-1].-Iu[2:N.x-1,2:N.y-1])/norm(Iu[2:N.x-1,2:N.y-1])
+ println("deviation: ",dev)
+
+ boundvalexp = [angle.(ψ[1:end-1,1]);angle.(ψ[end,1:end-1]);reverse(angle.(ψ[end,1:end-1]));reverse(angle.(ψ[1:end-1,1]))]
+
+ ## magnetic field ##
+ # magnetic induction
+ Bx = zeros(size(phase))
+ Bx[:, 2:end] = -hquer*N.x/(a.x*thick*qe) * diff(phase, dims=2)
+ Bx[:, 1] = Bx[:, 2]
+ By = zeros(size(phase))
+ By[2:end, :] = hquer * N.y / (a.y * thick * qe) * diff(phase, dims=1)
+ By[1, :] = Bx[2, :]
+
+ ## save results ##
+ if save
+ npzwrite("If_TIE.npy", If)
+ npzwrite("phase_TIE.npy", phase[2:end-1,2:end-1])
+ npzwrite("B_TIE.npy", Bx, By)
+ npzwrite("boundvalexp.npy", boundvalexp)
+ end
+end
+
+end # module
\ No newline at end of file
--
GitLab