From b24072b3befe6a6a0194ed72e7a4602caa3ef997 Mon Sep 17 00:00:00 2001 From: Carlos Castillo Passi Date: Wed, 26 Apr 2023 14:33:43 +0100 Subject: [PATCH] Fix for #162 #85 --- KomaMRICore/src/datatypes/Sequence.jl | 93 ++++++++++++++++++- .../src/simulation/TimeStepCalculation.jl | 22 +++-- 2 files changed, 104 insertions(+), 11 deletions(-) diff --git a/KomaMRICore/src/datatypes/Sequence.jl b/KomaMRICore/src/datatypes/Sequence.jl index f11b2f6ab..9ac143645 100644 --- a/KomaMRICore/src/datatypes/Sequence.jl +++ b/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 \ No newline at end of file diff --git a/KomaMRICore/src/simulation/TimeStepCalculation.jl b/KomaMRICore/src/simulation/TimeStepCalculation.jl index a902ac64b..fcf89b445 100644 --- a/KomaMRICore/src/simulation/TimeStepCalculation.jl +++ b/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,6 +240,7 @@ 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) @@ -244,8 +248,8 @@ function get_breaks_in_RF_key_points(seq::Sequence, t) 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