Merge branch 'koma_v074' of github.com:cncastillo/KomaMRI.jl into gui

KomaMRICore/src/datatypes/Sequence.jl

@@ -331,9 +331,9 @@ Generates a trapezoidal waveform vector.
 			B = A
 		end
 		#Trapezoidal waveform
-		aux = (ζ1    .<= t .- delay .< ζ1+ΔT) .* B
+		aux = (ζ1    .< t .- delay .< ζ1+ΔT) .* B
 		if ζ1 != 0
-			aux .+= (0     .<= t .- delay .< ζ1) .* A[1] .* (t .- delay)./ζ1
+			aux .+= (0     .< t .- delay .<= ζ1) .* A[1] .* (t .- delay)./ζ1
 		end
 		if ζ2 !=0
 			aux .+= (ζ1+ΔT .<= t .- delay .< ζ1+ΔT+ζ2) .* A[end] .* (1 .- (t.-delay.-ζ1.-ΔT)./ζ2)
@@ -813,4 +813,93 @@ get_M2(seq::Sequence; Δt=1) = begin
 	M2_adc = [M2x_adc M2y_adc M2z_adc]
 	#Final
 	M2, M2_adc
+end
+
+"""
+	SR, SR_adc = get_slew_rate(seq::Sequence; Δt=1)
+
+Outputs the designed slew rate of the Sequence `seq`.
+
+# Arguments
+- `seq`: (`::Sequence`) Sequence struct
+- `Δt`: (`::Real`, `=1`, `[s]`) nominal delta time separation between two time samples
+    for ADC acquisition and Gradients
+
+# Returns
+- `SR`: (`3-column ::Matrix{Float64}`) Slew rate
+- `SR_adc`: (`3-column ::Matrix{Float64}`) Slew rate sampled at ADC points
+"""
+get_slew_rate(seq::Sequence; Δt=1) = begin
+	t, Δt = get_uniform_times(seq, Δt)
+	Gx, Gy, Gz = get_grads(seq, t)
+	G = Dict(1=>[Gx; 0], 2=>[Gy; 0], 3=>[Gz; 0])
+	#kspace
+	Nt = length(t)
+	m2 = zeros(Nt,3)
+	#get_RF_center_breaks
+	idx_rf, rf_type = get_RF_types(seq, t)
+	parts = kfoldperm(Nt,1,type="ordered", breaks=idx_rf)
+	for i = 1:3
+		m2f = 0
+		for (rf, p) in enumerate(parts)
+			m2[p,i] = (G[i][p[1]+1:p[end]+1] .- G[i][p[1]:p[end]])[:] ./ Δt[p]
+		end
+	end
+	M2 = m2 #[m^-1]
+	#Interp, as Gradients are generally piece-wise linear, the integral is piece-wise cubic
+	#Nevertheless, the integral is sampled at the ADC times so a linear interp is sufficient
+	#TODO: check if this interpolation is necessary
+	ts = t .+ Δt
+	t_adc =  get_adc_sampling_times(seq)
+	M2x_adc = linear_interpolation(ts,M2[:,1],extrapolation_bc=0)(t_adc)
+	M2y_adc = linear_interpolation(ts,M2[:,2],extrapolation_bc=0)(t_adc)
+	M2z_adc = linear_interpolation(ts,M2[:,3],extrapolation_bc=0)(t_adc)
+	M2_adc = [M2x_adc M2y_adc M2z_adc]
+	#Final
+	M2, M2_adc
+end
+
+"""
+    EC, EC_adc = get_eddy_currents(seq::Sequence; Δt=1)
+
+Outputs the designed eddy currents of the Sequence `seq`.
+
+# Arguments
+- `seq`: (`::Sequence`) Sequence struct
+- `Δt`: (`::Real`, `=1`, `[s]`) nominal delta time separation between two time samples
+    for ADC acquisition and Gradients
+
+# Returns
+- `EC`: (`3-column ::Matrix{Float64}`) Eddy currents
+- `EC_adc`: (`3-column ::Matrix{Float64}`) Eddy currents sampled at ADC points
+"""
+get_eddy_currents(seq::Sequence; Δt=1, λ=80e-3) = begin
+	t, Δt = get_uniform_times(seq, Δt)
+	Gx, Gy, Gz = get_grads(seq, t)
+	G = Dict(1=>[Gx; 0], 2=>[Gy; 0], 3=>[Gz; 0])
+	#kspace
+	Nt = length(t)
+	m2 = zeros(Nt,3)
+	#get_RF_center_breaks
+	idx_rf, rf_type = get_RF_types(seq, t)
+	parts = kfoldperm(Nt,1,type="ordered", breaks=idx_rf)
+	for i = 1:3
+		m2f = 0
+		for (rf, p) in enumerate(parts)
+			m2[p,i] = (G[i][p[1]+1:p[end]+1] .- G[i][p[1]:p[end]])[:] ./ Δt[p]
+		end
+	end
+	ec(t, λ) = exp.(-t ./ λ) .* (t .>= 0)
+	M2 = [sum( m2[:, j] .* ec(t[i] .- t, λ) .* Δt ) for i = 1:Nt, j = 1:3] #[m^-1]
+	#Interp, as Gradients are generally piece-wise linear, the integral is piece-wise cubic
+	#Nevertheless, the integral is sampled at the ADC times so a linear interp is sufficient
+	#TODO: check if this interpolation is necessary
+	ts = t .+ Δt
+	t_adc =  get_adc_sampling_times(seq)
+	M2x_adc = linear_interpolation(ts,M2[:,1],extrapolation_bc=0)(t_adc)
+	M2y_adc = linear_interpolation(ts,M2[:,2],extrapolation_bc=0)(t_adc)
+	M2z_adc = linear_interpolation(ts,M2[:,3],extrapolation_bc=0)(t_adc)
+	M2_adc = [M2x_adc M2y_adc M2z_adc]
+	#Final
+	M2*1_000, M2_adc*1_000
 end

KomaMRICore/src/simulation/Bloch/BlochDictSimulationMethod.jl

@@ -3,6 +3,10 @@ Base.@kwdef struct BlochDict <: SimulationMethod
 end
 
 export BlochDict
+Base.show(io::IO, s::BlochDict) = begin
+	print(io, "BlochDict(save_Mz=$(s.save_Mz))")
+end
+
 
 output_Ndim(sim_method::BlochDict) = 3
 
@@ -80,8 +84,8 @@ It gives rise to a rotation of `M0` with an angle given by the efective magnetic
     a part of the complete Mag vector and it's a part of the initial state for the next
     precession simulation step)
 """
-function run_spin_excitation!(p::Phantom{T}, seq::DiscreteSequence{T}, 
+function run_spin_excitation!(p::Phantom{T}, seq::DiscreteSequence{T}, sig::AbstractArray{Complex{T}},
     M::Mag{T}, sim_method::BlochDict) where {T<:Real}
-    run_spin_excitation!(p, seq, M, Bloch()) #The same as Bloch
+    run_spin_excitation!(p, seq, sig, M, Bloch()) #The same as Bloch
     return nothing
 end

KomaMRICore/src/simulation/Bloch/BlochSimulationMethod.jl

@@ -76,7 +76,7 @@ It gives rise to a rotation of `M0` with an angle given by the efective magnetic
     a part of the complete Mag vector and it's a part of the initial state for the next
     precession simulation step)
 """
-function run_spin_excitation!(p::Phantom{T}, seq::DiscreteSequence{T}, 
+function run_spin_excitation!(p::Phantom{T}, seq::DiscreteSequence{T}, sig::AbstractArray{Complex{T}},
     M::Mag{T}, sim_method::Bloch) where {T<:Real}
     #Simulation
     for s ∈ seq #This iterates over seq, "s = seq[i,:]"
@@ -96,5 +96,7 @@ function run_spin_excitation!(p::Phantom{T}, seq::DiscreteSequence{T},
         M.xy .= M.xy .* exp.(-s.Δt ./ p.T2)
         M.z  .= M.z  .* exp.(-s.Δt ./ p.T1) .+ p.ρ .* (1 .- exp.(-s.Δt ./ p.T1))
     end
+    #Acquired signal
+    #sig .= -0.1im #<-- This was to test if an ADC point was inside an RF block
     return nothing
 end

KomaMRICore/src/simulation/SimulatorCore.jl

@@ -58,14 +58,15 @@ different number threads to excecute the process.
 - `M0`: (`::Vector{Mag}`) final state of the Mag vector after a rotation (or the initial
     state for the next precession simulation step)
 """
-function run_spin_excitation_parallel!(obj::Phantom{T}, seq::DiscreteSequence{T},
+function run_spin_excitation_parallel!(obj::Phantom{T}, seq::DiscreteSequence{T}, sig::AbstractArray{Complex{T}},
     Xt::SpinStateRepresentation{T}, sim_method::SimulationMethod;
     Nthreads=Threads.nthreads()) where {T<:Real}
 
     parts = kfoldperm(length(obj), Nthreads; type="ordered")
+    dims = [Colon() for i=1:output_Ndim(sim_method)] # :,:,:,... Ndim times
 
-    Threads.@threads for p ∈ parts
-        run_spin_excitation!(@view(obj[p]), seq, @view(Xt[p]), sim_method)
+    ThreadsX.foreach(enumerate(parts)) do (i, p)
+        run_spin_excitation!(@view(obj[p]), seq, @view(sig[dims...,i]), @view(Xt[p]), sim_method)
     end
 
     return nothing
@@ -112,12 +113,12 @@ function run_sim_time_iter!(obj::Phantom, seq::DiscreteSequence, sig::AbstractAr
         dims = [Colon() for i=1:output_Ndim(sim_method)] # :,:,:,... Ndim times
         # Simulation wrappers
         if excitation_bool
-            run_spin_excitation_parallel!(obj, seq_block, Xt, sim_method; Nthreads)
+            run_spin_excitation_parallel!(obj, seq_block, @view(sig[acq_samples, dims...]), Xt, sim_method; Nthreads)
             rfs += 1
         else
             run_spin_precession_parallel!(obj, seq_block, @view(sig[acq_samples, dims...]), Xt, sim_method; Nthreads)
-            samples += Nadc
         end
+        samples += Nadc
         #Update progress
         next!(progress_bar, showvalues=[(:simulated_blocks, block), (:rf_blocks, rfs), (:acq_samples, samples-1)])
         update_blink_window_progress!(w, block, Nblocks)

KomaMRICore/src/simulation/TimeStepCalculation.jl

@@ -1,3 +1,4 @@
+const EPS = eps(1.0)
 """
     array_of_ranges = kfoldperm(N, k; type="random", breaks=[])
 
@@ -178,23 +179,19 @@ This function returns non-uniform time points that are relevant in the sequence
 """
 function get_variable_times(seq; dt=1e-3, dt_rf=1e-5)
 	t = Float64[]
+	ϵ = EPS #Smallest Float64
 	ΔT = durs(seq) #Duration of sequence block
 	T0 = cumsum([0; ΔT[:]]) #Start time of each block
 	for i = 1:length(seq)
 		s = seq[i] #Current sequence block
 		t0 = T0[i]
-		if is_ADC_on(s)
-			ts = get_adc_sampling_times(s) .+ t0
-			taux = points_from_key_times(ts; dt) # ADC sampling
-			append!(t, taux)
-		end
 		if is_RF_on(s)
 			y = s.RF[1]
 			delay, T = y.delay, y.T
 			t1 = t0 + delay
 			t2 = t1 + sum(T)
 			rf0 = t0 + get_RF_center(y) #get_RF_center includes delays
-			taux = points_from_key_times([t1,rf0,t2]; dt=dt_rf) # Arbitrary RF
+			taux = points_from_key_times([t1,t1+ϵ,rf0,t2-ϵ,t2]; dt=dt_rf) # Arbitrary RF. Points (t1+ϵ, t2-ϵ) added to fix bug with ADCs
 			append!(t, taux)
 		end
 		if is_GR_on(s)
@@ -204,13 +201,19 @@ function get_variable_times(seq; dt=1e-3, dt_rf=1e-5)
 			if is_Gz_on(s) append!(active_gradients, s.GR.z) end
 			for y = active_gradients
 				ts = get_theo_t(y) .+ t0
-				taux = points_from_key_times(ts; dt)
+				taux = points_from_key_times([ts[1]+ϵ; ts; ts[end]-ϵ]; dt) #The ±ϵ fixes #
 				append!(t, taux)
 			end
 		end
 	end
+	#Adding ADC samples, and removing repeated points
 	tadc = get_adc_sampling_times(seq)
 	t = sort(unique([t; tadc])) #Removing repeated points
+	#Fixes a problem with ADC at the start and end of the seq
+	t0 = t[1]   - ϵ
+	tf = t[end] + ϵ
+	t = [t0; t; tf]
+	#Final time points
 	Δt = t[2:end] .- t[1:end-1]
 	t = t[1:end-1]
 	t, Δt
@@ -237,15 +240,16 @@ Thus, it is possible to split the simulation into excitation and preccesion comp
 function get_breaks_in_RF_key_points(seq::Sequence, t)
 	T0 = cumsum([0; durs(seq)[:]])
 	# Identifying RF key points
+	ϵ = EPS #Smallest Float64
 	key_points = Float64[]
 	key_idxs = Int[]
 	for (i, s) = enumerate(seq)
 		if is_RF_on(s)
 			t0 = T0[i] + s.RF.delay[1]	#start of RF waverform
 			tf = T0[i] + s.RF.dur[1]	#end of RF waveform
 			append!(key_points, [t0; tf])
-			idx0 = argmin(abs.(t.-t0))
-			idxf = argmin(abs.(t.-tf))
+			idx0 = argmin(abs.(t.-(t0+ϵ)))
+			idxf = argmin(abs.(t.-(tf-ϵ)))
 			append!(key_idxs,   [idx0; idxf])
 		end
 	end

KomaMRIPlots/src/KomaMRIPlots.jl

@@ -8,6 +8,6 @@ include("ui/DisplayFunctions.jl")
 using Reexport
 @reexport using PlotlyJS: savefig
 
-export plot_seq, plot_M0, plot_M1, plot_M2, plot_kspace, plot_phantom_map, plot_signal, plot_image, plot_dict
+export plot_seq, plot_M0, plot_M1, plot_M2, plot_eddy_currents, plot_kspace, plot_phantom_map, plot_signal, plot_image, plot_dict
 
 end

KomaMRIPlots/src/ui/DisplayFunctions.jl

@@ -91,7 +91,7 @@ julia> seq = read_seq(seq_file)
 julia> plot_seq(seq)
 ```
 """
-plot_seq(seq::Sequence; width=nothing, height=nothing, slider=true, show_seq_blocks=false, show_sim_blocks=false, Nblocks=0, darkmode=false, max_rf_samples=100, range=[]) = begin
+plot_seq(seq::Sequence; width=nothing, height=nothing, slider=true, show_seq_blocks=false, show_sim_blocks=false, Nblocks=0, darkmode=false, max_rf_samples=100, range=[], title="") = begin
 	bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode)
 	idx = ["Gx" "Gy" "Gz"]
 	N = length(seq)
@@ -129,7 +129,7 @@ plot_seq(seq::Sequence; width=nothing, height=nothing, slider=true, show_seq_blo
 		t, _ = KomaMRICore.get_uniform_times(seq, 1e-3)
 		breaks = KomaMRICore.get_breaks_in_RF_key_points(seq,t)
 		Nt = length(t)
-		if Nblocks == 0
+		if Nblocks == 0 #TODO: This should change to a call to a function that generates the default parameters for the simulation 
 			Nblocks = 20
 		end
 		parts = KomaMRICore.kfoldperm(Nt,Nblocks;type="ordered",breaks)
@@ -143,7 +143,7 @@ plot_seq(seq::Sequence; width=nothing, height=nothing, slider=true, show_seq_blo
 			) for i = 1:length(t_sim_parts)]
 		append!(shapes, aux)
 	end
-	l = PlotlyJS.Layout(; hovermode="closest",
+	l = PlotlyJS.Layout(;title=title, hovermode="closest",
 			xaxis_title="",
 			modebar=attr(orientation="h",yanchor="bottom",xanchor="right",y=1,x=0,bgcolor=bgcolor,color=text_color,activecolor=plot_bgcolor),
 			legend=attr(orientation="h",yanchor="bottom",xanchor="left",y=1,x=0),
@@ -515,7 +515,7 @@ end
 
 
 """
-    p = plot_M1(seq; height=nothing, width=nothing, slider=true, darkmode=false)
+    p = plot_M2(seq; height=nothing, width=nothing, slider=true, darkmode=false)
 
 Plots the second order moment (M2) of a Sequence `seq`.
 
@@ -601,6 +601,94 @@ function plot_M2(seq; height=nothing, width=nothing, slider=true, darkmode=false
 	PlotlyJS.plot(p, l; config)
 end
 
+
+"""
+    p = plot_eddy_currents(seq; height=nothing, width=nothing, slider=true, darkmode=false)
+
+Plots the eddy currents of a Sequence `seq`.
+
+# Arguments
+- `seq`: (`::Sequence`) Sequence
+
+# Keywords
+- `height`: (`::Int64`, `=nothing`) height of the plot
+- `width`: (`::Int64`, `=nothing`) width of the plot
+- `slider`: (`::Bool`, `=true`) boolean to display a slider
+- `darkmode`: (`::Bool`, `=false`) boolean to define colors for darkmode
+
+# Returns
+- `p`: (`::PlotlyJS.SyncPlot`) plot of the eddy currents of the sequence struct `seq`
+
+# Examples
+```julia-repl
+julia> seq_file = joinpath(dirname(pathof(KomaMRI)), "../examples/1.sequences/spiral.seq")
+
+julia> seq = read_seq(seq_file)
+
+julia> plot_eddy_currents(seq)
+```
+"""
+function plot_eddy_currents(seq, λ; height=nothing, width=nothing, slider=true, darkmode=false, range=[])
+	#Times
+	dt = 1
+	t, Δt = KomaMRICore.get_uniform_times(seq, dt)
+	#kx,ky
+	ts = t .+ Δt
+	rf_idx, rf_type = KomaMRICore.get_RF_types(seq, t)
+	k, _ =  KomaMRICore.get_eddy_currents(seq; Δt=dt, λ)
+
+	#plots k(t)
+	bgcolor, text_color, plot_bgcolor, grid_color, sep_color = theme_chooser(darkmode)
+	l = PlotlyJS.Layout(;yaxis_title="EC", hovermode="closest",
+			xaxis_title="",
+			modebar=attr(orientation="h",yanchor="bottom",xanchor="right",y=1,x=0,bgcolor=bgcolor,color=text_color,activecolor=plot_bgcolor),
+			legend=attr(orientation="h",yanchor="bottom",xanchor="left",y=1,x=0),
+			plot_bgcolor=plot_bgcolor,
+			paper_bgcolor=bgcolor,
+			xaxis_gridcolor=grid_color,
+			yaxis_gridcolor=grid_color,
+			xaxis_zerolinecolor=grid_color,
+			yaxis_zerolinecolor=grid_color,
+			font_color=text_color,
+			yaxis_fixedrange = false,
+			xaxis=attr(
+				ticksuffix=" ms",
+				range=range[:],
+				rangeslider=attr(visible=slider),
+				rangeselector=attr(
+					buttons=[
+						attr(count=1,
+						label="1m",
+						step=10,
+						stepmode="backward"),
+						attr(step="all")
+						]
+					),
+				),
+			margin=attr(t=0,l=0,r=0)
+			)
+    if height !== nothing
+        l.height = height
+    end
+    if width !== nothing
+        l.width = width
+    end
+	plotter = PlotlyJS.scatter #using scattergl speeds up the plotting but does not show the sequence in the slider below
+	p = [plotter() for j=1:4]
+	p[1] = plotter(x=ts*1e3, y=k[:,1], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms³)", name="M2x")
+	p[2] = plotter(x=ts*1e3, y=k[:,2], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms³)", name="M2y")
+	p[3] = plotter(x=ts*1e3, y=k[:,3], hovertemplate="(%{x:.4f} ms, %{y:.2f} mT/m⋅ms³)", name="M2z")
+	# p[4] = plotter(x=t[rf_idx]*1e3,y=rf_type,name="RFs",marker=attr(symbol="cross",size=8,color="orange"),mode="markers")
+	config = PlotConfig(
+		displaylogo=false,
+		toImageButtonOptions=attr(
+			format="svg", # one of png, svg, jpeg, webp
+		).fields,
+		modeBarButtonsToRemove=["zoom", "select", "select2d","lasso2d", "autoScale", "resetScale2d", "pan", "tableRotation", "resetCameraLastSave", "zoomIn", "zoomOut"]
+	)
+	PlotlyJS.plot(p, l; config)
+end
+
 """
     p = plot_phantom_map(ph, key; t0=0, height=600, width=nothing, darkmode=false)
 

README.md

@@ -89,6 +89,7 @@ To install just write the following in the Julia REPL:
 KomaMRI.jl comes with a handy GUI that contains a brain phantom with an EPI sequence. To open it use:
 
 ```julia
+using KomaMRI
 KomaUI()
 ```
 Press the button that says "Simulate!" to do your first simulation :). Then, a notification emerge telling you that the simulation was successful. In this notification, you can either select to (1) see the Raw signal or (2) to procced with the reconstruction.