created BlochVoxelDictSimulation simulation method

KomaMRICore/src/simulation/Bloch/BlochVoxelDictSimulation.jl

@@ -0,0 +1,50 @@
+import KomaMRICore: SimulationMethod, initialize_spins_state, mul!, output_Ndim,
+    run_spin_excitation!, run_spin_precession!, sim_output_dim
+
+
+Base.@kwdef struct BlochVoxelDict <: SimulationMethod 
+    N_spins_per_voxel::Int=1
+end
+
+
+function sim_output_dim(obj::Phantom{T}, seq::Sequence, sys::Scanner, sim_method::BlochVoxelDict) where {T<:Real}
+    N_voxels = length(obj) ÷ sim_method.N_spins_per_voxel
+    return (sum(seq.ADC.N), N_voxels)
+end
+
+
+function run_spin_precession!(p::Phantom{T}, seq::DiscreteSequence{T}, sig::AbstractArray{Complex{T}}, 
+    M::Mag{T}, sim_method::BlochVoxelDict) where {T<:Real}
+    #Simulation
+    #Motion
+    xt = p.x .+ p.ux(p.x, p.y, p.z, seq.t')
+    yt = p.y .+ p.uy(p.x, p.y, p.z, seq.t')
+    zt = p.z .+ p.uz(p.x, p.y, p.z, seq.t')
+    #Effective field
+    Bz = xt .* seq.Gx' .+ yt .* seq.Gy' .+ zt .* seq.Gz' .+ p.Δw / T(2π * γ)
+    #Rotation
+    if is_ADC_on(seq)
+        ϕ = T(-2π * γ) .* cumtrapz(seq.Δt', Bz)
+    else
+        ϕ = T(-2π * γ) .* trapz(seq.Δt', Bz)
+    end
+    #Mxy precession and relaxation, and Mz relaxation
+    tp = cumsum(seq.Δt) # t' = t - t0
+    dur = sum(seq.Δt)   # Total length, used for signal relaxation
+    Mxy = [M.xy M.xy .* exp.(1im .* ϕ .- tp' ./ p.T2)] #This assumes Δw and T2 are constant in time
+    M.xy .= Mxy[:, end]
+    M.z .= M.z .* exp.(-dur ./ p.T1) .+ p.ρ .* (1 .- exp.(-dur ./ p.T1))
+
+    #Acquired signal
+    # get indexes of the voxels being processed by the present thread
+    idx = floor.(Int, p.Dλ1)
+    # for each voxel in the voxel-array (list of phantoms) assign the corresponding magnetization value (sum over the spins within a voxel) and normalize
+    N_voxels = size(sig, 2)
+    #=     println(sum(Mxy[findall(idx .== 1), findall(seq.ADC)];dims=1)/ sim_method.N_spins_per_voxel )
+    sleep(3) =#
+    for voxel in 1:N_voxels
+        sig[:, voxel] += transpose(sum(Mxy[findall(idx .== voxel), findall(seq.ADC)];dims=1)) / sim_method.N_spins_per_voxel 
+    end
+
+    return nothing
+end

KomaMRIFiles/src/Phantom/T1T2_voxel_phantom.jl

@@ -0,0 +1,33 @@
+using KomaMRI
+
+
+B0 = 0.55         # 0.55 T
+Niso = 1          # number of spins in x and y directions 
+Niso_z = 101      # number of spins in z direction
+Δx_voxel = 2e-3   # x voxel dim 
+Δy_voxel = 2e-3   # y voxel dim
+Δz_voxel = 10e-3  # z voxel dim
+
+NisoTotal = Niso*Niso*Niso_z
+
+if Niso == 1
+    dx = 0 
+    dy = 0
+else
+    dx = Array(range(-Δx_voxel/2, Δx_voxel/2, Niso))
+    dy = Array(range(-Δy_voxel/2, Δy_voxel/2, Niso))
+end
+
+dz = Array(range(-Δz_voxel/2, Δz_voxel/2, Niso_z))
+coords = vec(collect(Iterators.product(dx, dy, dz)))
+
+function voxel_phantom(T1v, T2v, idx)
+    obj = Phantom{Float64}(x=[x[1] for x in coords], 
+                           y=[x[2] for x in coords],
+                           z=[x[3] for x in coords],
+                           T1=T1v*ones(Niso*Niso*Niso_z), 
+                           T2=T2v*ones(Niso*Niso*Niso_z),
+                           Dλ1=idx*ones(Niso*Niso*Niso_z))
+    return obj
+end
+

KomaMRIFiles/src/Phantom/t1mes_lowfield_phantom.jl

@@ -0,0 +1,86 @@
+using KomaMRIPlots, KomaMRICore, Distributions
+
+X = collect(1:128) 
+Y = collect(1:128)
+
+positions = [(32, 32), (64, 32), (96, 32), (32, 64), (64, 64), (96, 64), (32, 96), (64, 96), (96, 96)]
+
+r = 11 # 22mm diameter per vial
+M = zeros(UInt32, 128, 128, 9) 
+for i in 1:9
+    pos_x, pos_y = positions[i]
+    M[:,:,i] .= i.*((X .- pos_x).^2 .+ (Y' .- pos_y).^2 .<= r.^2)
+end
+
+M = sum(M, dims=3)
+M = convert(Array{UInt32, 3}, M)
+class = Matrix{UInt32}(reshape(M, 128, 128))'
+
+B0 = 0.55         # 0.55 T
+Niso = 200        # 200 isochromats
+
+# Define spin position arrays
+Δx = (122/117)*.1e-2        # 10mm
+Δz_voxel = 10e-3            # 10 mm
+M, N = size(class)          # Number of spins in x and y
+FOVx = (M-1)*Δx             # Field of view in x
+FOVy = (N-1)*Δx             # Field of view in y
+x = -FOVx/2:Δx:FOVx/2       # x spin coordinates vector
+y = -FOVy/2:Δx:FOVy/2       # y spin coordinates vector
+x, y = x .+ y'*0, x*0 .+ y' # x and y grid points
+
+
+# Define the proton density array
+ρ = (class.== 0)*1 .+    # Matrix
+    (class.== 1)*1 .+    # vial 1
+    (class.== 2)*1 .+    # vial 2
+    (class.== 3)*1 .+    # vial 3 
+    (class.== 4)*1 .+    # vial 4
+    (class.== 5)*1 .+    # vial 5
+    (class.== 6)*1 .+    # vial 6
+    (class.== 7)*1 .+    # vial 7
+    (class.== 8)*1 .+    # vial 8
+    (class.== 9)*1       # vial 9
+
+# Define the T1 decay array
+import Random 
+Random.seed!(1234)
+T1 = (class.==0)*810 .+ (class.==0).*Random.rand(Uniform(-16,16),M,N) .+  # Matrix
+     (class.==1)*424 .+ (class.==1).*Random.rand(Uniform(-4,4),M,N) .+    # vial 1
+     (class.==2)*993 .+ (class.==2).*Random.rand(Uniform(-4,4),M,N) .+    # vial 2
+     (class.==3)*450 .+ (class.==3).*Random.rand(Uniform(-5,5),M,N) .+    # vial 3
+     (class.==4)*543 .+ (class.==4).*Random.rand(Uniform(-3,3),M,N) .+    # vial 4
+     (class.==5)*1185 .+ (class.==5).*Random.rand(Uniform(-7,7),M,N) .+   # vial 5
+     (class.==6)*1460 .+ (class.==6).*Random.rand(Uniform(-15,15),M,N) .+ # vial 6
+     (class.==7)*299 .+ (class.==7).*Random.rand(Uniform(-3,3),M,N) .+    # vial 7
+     (class.==8)*751 .+ (class.==8).*Random.rand(Uniform(-4,4),M,N) .+    # vial 8
+     (class.==9)*264 .+ (class.==9).*Random.rand(Uniform(-3,3),M,N)       # vial 9
+
+# Define the T2 decay array
+T2 = (class.==0)*125 .+ (class.==0).*Random.rand(Uniform(-20,20),M,N) .+  # Matrix
+     (class.==1)*47 .+ (class.==1).*Random.rand(Uniform(-1,1),M,N) .+     # vial 1
+     (class.==2)*53 .+ (class.==2).*Random.rand(Uniform(-1.5,1.5),M,N) .+ # vial 2
+     (class.==3)*196 .+ (class.==3).*Random.rand(Uniform(-5,5),M,N) .+    # vial 3
+     (class.==4)*49 .+ (class.==4).*Random.rand(Uniform(-1,1),M,N) .+     # vial 4
+     (class.==5)*53 .+ (class.==5).*Random.rand(Uniform(-1,1),M,N) .+     # vial 5
+     (class.==6)*238 .+ (class.==6).*Random.rand(Uniform(-5,5),M,N) .+    # vial 6
+     (class.==7)*47 .+ (class.==7).*Random.rand(Uniform(-1,1),M,N) .+     # vial 7
+     (class.==8)*53 .+ (class.==8).*Random.rand(Uniform(-1,1),M,N) .+     # vial 8
+     (class.==9)*174 .+ (class.==9).*Random.rand(Uniform(-3,3),M,N)       # vial 9
+
+# Time scaling factor
+T1 = T1*1e-3
+T2 = T2*1e-3
+
+# Define the phantom object
+obj = Phantom{Float64}(
+    name = "custom_T1MES_lowfield",
+	x = x[ρ.!=0],
+	y = y[ρ.!=0],
+	z = 0*x[ρ.!=0],
+	ρ = ρ[ρ.!=0],
+	T1 = T1[ρ.!=0],
+	T2 = T2[ρ.!=0],
+)
+
+

examples/dict_sim_example.jl

@@ -0,0 +1,93 @@
+using KomaMRI, PlotlyJS, MAT
+
+# 1. DEFINE MAIN PARAMETERS
+# 1.1 Define lists for tissue parameters - this is just an example of possible values
+T1list = Array(vcat(50:50:500, 505:5:1000, 1050:50:1500, 1600:100:2000))*1e-3
+T2list = Array(vcat(4:1:80, 85:5:200, 220:20:400))*1e-3
+
+# 1.2 Compute combinations (T1-T2 pairs)
+params_combs = vec(collect(Iterators.product(T1list, T2list, B1list, T1rholist, Dlist)))
+# get ordered lists for T1-T2 values (for consistency)
+T1values = [x[1] for x in params_combs]
+T2values = [x[2] for x in params_combs]
+Ncombinations = length(params_combs)
+
+
+# 2. READ/LOAD SEQUENCE
+seq_path = "path/to/your_seq.seq"
+seq = read_seq(seq_path)
+
+# 2.1 Add heart rate variability - if needed
+# read .mat file with RRs vector
+# RRs_struct is a MATLAB struct with a field subject_RRs which contains an array with the measured RRs
+# this is just a way to do it, it may be different depending on your application
+RRs_struct = matread("MRF_simulations/RRs/HV8_RRs.mat")
+desired_RRs = RRs_struct["subject_RRs"][:]*1e-3
+
+# 2.2 Get triggered blocks from sequence
+trigg_idx = findall(seq.DEF["extension"] .!= 0)
+t0 = get_block_start_times(seq)[1:end-1][trigg_idx] # The 1:end-1 is because start time are [0, t0, t1, .., tn]
+current_RRs = diff(t0)
+
+# 2.3 Modify RRs to the actual value
+seq.DUR[trigg_idx[1:end-1]] .+= (desired_RRs .- current_RRs)
+t0 = get_block_start_times(seq)[1:end-1][trigg_idx] 
+new_RRs = diff(t0)
+
+## plot the sequence - if you want
+using KomaMRIPlots.PlotlyJS
+p = plot_seq(seq)
+for i in eachindex(t0)
+    add_trace!(p, KomaMRIPlots.scatter(x=[t0[i], t0[i]]*1e3, y=[-30, 30],
+        marker=attr(color="red"), legendgroup="RRs", showlegend=i==1, name="RRs")) # we mark the RR timestamps
+end
+p
+
+# 2.4 Only sample the center of the readout - you may need to modify this for sequences with echoes
+# you could also add more samples but computation time increases significantly
+for (i, adc) in enumerate(seq.ADC[is_ADC_on.(seq)])
+    adc.N = 1
+end
+
+
+# 3. SIMULATION
+# 3.1 Define main simulation parameters
+include("path/to/T1T2_voxel_phantom.jl")
+include("path/to/BlochVoxelDictSimulation.jl")
+sim_params = KomaMRICore.default_sim_params()
+sim_params["sim_method"] =  BlochVoxelDict(N_spins_per_voxel=NisoTotal) # N_spins_per_voxel depends on the spins defined in the voxelphantom 
+sim_params["Nthreads"] = 8        # define number of threads - depends on your machine
+sim_params["return_type"] = "mat" # format for the output signal
+sim_params["gpu"] = false         # GPU not supported yet
+
+# system struct instance
+sys = Scanner()
+
+# 3.2 Generate list of phantoms - each voxel (phantom) has a unique T1-T2 combination and a unique ID (n)
+phantom_list = []
+for n in range(1,Ncombinations)
+    push!(phantom_list, voxel_phantom(T1values[n],T2values[n], n)) 
+end
+
+# 3.3 sum list of phantoms and run the simulation!
+phantom = sum(phantom_list)
+raw = simulate(phantom, seq, sys; sim_params)
+
+
+# 4. PLOT STUFF
+# maybe you want to plot the results, e.g. plot the magnitude of the first 5 entries of the dictionary
+sig = abs.(raw[:,1:5,1])
+plot(sig)
+
+
+# 5. SAVE RESULTS
+# reshape combinations lists as an array
+lists_temp = [[1e3*k for k in x] for x in params_combs] # in ms
+lists_combinations = mapreduce(permutedims, vcat, lists_temp)
+
+# save MATLAB struct named KomaDict with fields
+#                          KomaDict.raw
+#                          KomaDict.combinations
+# you can also add formatting for a more descriptive name
+save_folder = "path/to/save/"
+matwrite(join([save_folder, "KomaDict.mat"]), Dict("GenDict" => (Combinations = lists_combinations, dictOn = raw)))