Upload New File

src/TIE.jl

+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