diff --git a/julia/pseudoinverse.jl b/julia/pseudoinverse.jl index 3a6b683..743ce1a 100644 --- a/julia/pseudoinverse.jl +++ b/julia/pseudoinverse.jl @@ -1,40 +1,27 @@ -@doc """ -This type stores the singular value decomposition V⁺ΣU of a matrix A as well -as the inverse singular values, which are neccessary to calculate its -pseudoinverse. -"""-> -type SVD - U::Matrix - Σ::Vector - V::Matrix - D::Vector -end - -SVD(U::Matrix,Σ::Vector,V::Matrix) = SVD(U,Σ,V,1./Σ) - @doc """ This algorithm solves the Thikonov regularized least squares Problem argminₓ(‖Ax-b‖² + λ‖b‖²) using the singular value decomposition of A. # Arguments -* `SVD::SVD`: Singular value decomposition of A +* `U, Σ, V`: Singular value decomposition of A * `b::Vector`: Measurement vector b * `lambd::Float64`: The regularization parameter, relative to the matrix trace * `enforceReal::Bool`: Enable projection of solution on real plane during iteration * `enforcePositive::Bool`: Enable projection of solution onto positive halfplane during iteration """ -> -function pseudoinverse{T}(S::SVD, b::Vector{T}, lambd, enforceReal, enforcePositive) +function pseudoinverse{T}(U::Matrix, Σ::Vector, V::Matrix, b::Vector{T}, lambd, enforceReal, enforcePositive) # perform regularization - for i=1:length(S.Σ) - σi = S.Σ[i] - S.D[i] = σi/(σi^2+lambd^2) + D = zeros(Σ) + for i=1:length(Σ) + σi = Σ[i] + D[i] = σi/(σi^2+lambd^2) end # calculate pseudoinverse - tmp = BLAS.gemv('C', one(T), S.U, b) - tmp .*= S.D - c = BLAS.gemv('N', one(T), S.V, tmp) + tmp = BLAS.gemv('C', one(T), U, b) + tmp .*= D + c = BLAS.gemv('N', one(T), V, tmp) # apply constraints if enforceReal && eltype(c) <: Complex diff --git a/julia/reco.jl b/julia/reco.jl index 3449029..6a4da06 100644 --- a/julia/reco.jl +++ b/julia/reco.jl @@ -40,8 +40,8 @@ u = vec(mean(u,2)) c = kaczmarzReg(S,u,1,1e6,false,true,true) # reconstruct using signular value decomposition -SSVD = SVD(svd(S.')...) -csvd = pseudoinverse(SSVD, u, 5e3, true, true) +U, Σ, V = svd(S.') +csvd = pseudoinverse(U, Σ, V, u, 5e3, true, true) # reshape into an image N = h5read(filenameSM, "/calibration/size") @@ -51,9 +51,9 @@ csvd = reshape(csvd,N[1],N[2]) # plot kaczmarz reconstruction figure() gray() -imshow(real(c)) +imshow(real(c), interpolation="None") # plot pseudoinverse reconstruction figure() gray() -imshow(real(csvd)) +imshow(real(csvd), interpolation="None") diff --git a/specification/MDF.pdf b/specification/MDF.pdf index a43e34e..29916c9 100644 Binary files a/specification/MDF.pdf and b/specification/MDF.pdf differ diff --git a/specification/MDF.tex b/specification/MDF.tex index 1f18197..db00ef5 100644 --- a/specification/MDF.tex +++ b/specification/MDF.tex @@ -64,16 +64,16 @@ \begin{document} \title{MDF: Magnetic Particle Imaging Data Format} -\newcommand{\version}{1.0.2} +\newcommand{\version}{1.0.3} \author{ -T.~Knopp$^{1,2}$, T.~Viereck$^{3}$, G.~Bringout$^{5}$, M.~Ahlborg$^{6}$, J.~Rahmer$^7$, M.~Hofmann$^{1,2}$ \\ \\ +T.~Knopp$^{1,2}$, T.~Viereck$^{3}$, G.~Bringout$^{4}$, M.~Ahlborg$^{5}$, J.~Rahmer$^6$, M.~Hofmann$^{1,2}$ \\ \\ $^1$Section for Biomedical Imaging, University Medical Center Hamburg-Eppendorf, Germany\\ $^2$Institute for Biomedical Imaging, Hamburg University of Technology, Germany\\ $^3$Institute of Electrical Measurement and Fundamental Electrical Engineering, TU Braunschweig, Germany\\ -$^5$Physikalisch-Technische Bundesanstalt, Berlin, Germany\\ -$^6$Institute of Medical Engineering, University of Lübeck, Germany\\ -$^7$Philips GmbH Innovative Technologies, Research Laboratories, Röntgenstraße 24-26, 22315 Hamburg, Germany +$^4$Physikalisch-Technische Bundesanstalt, Berlin, Germany\\ +$^5$Institute of Medical Engineering, University of Lübeck, Germany\\ +$^6$Philips GmbH Innovative Technologies, Research Laboratories, Röntgenstraße 24-26, 22315 Hamburg, Germany } %\address{