Merge pull request #30 from hofmannmartin/master

update to 1.0.3

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

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")

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{