Open katyabrui opened 8 months ago
@katyabrui
Great that you are using KomaMRI, and thanks for reporting the issue!
An immediate work around might be to use:
your_obj = brain_phantom3D(; ss=your_ss, start_end=[177, 183]);
if using KomaUI:
obj_ui[] = brain_phantom3D(; ss=your_ss, start_end=[177, 183]);
This is a 3D phantom that is very thin along Z, will take a bit longer to simulate, but less so than the default brain_phantom3D(). ss defaults to 4 for both 2D and 3D.
Hope this helps while the issue is investigated!
=====================
This unfortunately does not fix the issue.
@katyabrui
I am trying the seq now. Nice job getting it through github!
This is related to #308, but instead of the "echo dimension," the problem is in the "slice dimension." Right now, we are guessing how to reconstruct stuff (the dimensions Nx, Ny, Nz); a proper way would be to specify the label counters in the "[EXTENSION]" section of the pulseq file. I believe previously, there was a warning that told the user we were doing this. We should put it back (@beorostica). For now, for the echoes, we could add a seq.DEF["Necho"]
to reconstruct interleaved echoes.
In the meantime, you could specify the correct sequence dimensions in the REPL before reconstructing.
seq_ui[].DEF["Nx"] = 32
seq_ui[].DEF["Ny"] = 32
seq_ui[].DEF["Nz"] = 32
I can not test it right now, but I think it should work. For convenience, this also can be specified in the sequence file
[DEFINITIONS]
AdcRasterTime 6.25e-05
BlockDurationRaster 1e-05
GradientRasterTime 1e-05
RadiofrequencyRasterTime 1e-06
Nx 32
Ny 32
Nz 32
@cncastillo @katyabrui
This seems to be the issue, the current seq_ui[].DEF below has Nz=1:
julia> seq_ui[].DEF
Dict{String, Any} with 12 entries:
"signature" => "9adc63c8e1f8331d67b3fb1c9bdeac98"
"AdcRasterTime" => 6.25e-5
"GradientRasterTime" => 1.0e-5
"Nz" => 1
"Nx" => 32
"Ny" => 1024
"PulseqVersion" => 1004000
"BlockDurationRaster" => 1.0e-5
"extension" => [0, 0, 0, 0, 0, 0, 0, 0, 0, 0 … 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
"FileName" => "MPRAGE_brain_25cmfov.seq"
"RadiofrequencyRasterTime" => 1.0e-6
"additional_text" => "# Format of extension lists:\n# id type ref next_id\n# next_id of 0 terminates the list\n# Extension list is followed by extension specification…
The Nz, Nx, Ny parts of DEF are generated by KomaMRI(?) and are not in the current seq file:
[DEFINITIONS]
AdcRasterTime 6.25e-05
BlockDurationRaster 1e-05
GradientRasterTime 1e-05
RadiofrequencyRasterTime 1e-06
I saw a similar issue with:
https://github.com/pulseq/pulseq/blob/master/matlab/demoSeq/writeGradientEcho3D.m
but had not gotten around to resolving yet. Will try the above fix, thanks!
@cncastillo @katyabrui
Changing the Nx-z parameters in seq.DEF:
julia> seq.DEF
Dict{String, Any} with 12 entries:
"signature" => "9adc63c8e1f8331d67b3fb1c9bdeac98"
"AdcRasterTime" => 6.25e-5
"GradientRasterTime" => 1.0e-5
"Nz" => 32
"Nx" => 32
"Ny" => 32
"PulseqVersion" => 1004000
"BlockDurationRaster" => 1.0e-5
"extension" => [0, 0, 0, 0, 0, 0, 0, 0, 0, 0 … 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
"FileName" => "MPRAGE_brain_25cmfov.seq"
"RadiofrequencyRasterTime" => 1.0e-6
"additional_text" => "# Format of extension lists:\n# id type ref next_id\n# next_id of 0 terminates the list\n# Extension list is followed by extension specifications\n[EXTENSIONS]\n\n#…
And running the simulation, results in:
julia> raw = simulate(obj, seq, sys);
┌ Info: Running simulation in the GPU (Quadro RTX 8000)
│ koma_version = v"0.8.1"
│ sim_method = Bloch()
│ spins = 78126
│ time_points = 81897
└ adc_points = 32768
Progress: 100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:01:04
simulated_blocks: 2264
rf_blocks: 1124
acq_samples: 32768
65.210778 seconds (65.12 M allocations: 3.968 GiB, 2.15% gc time)
which seems to be simulating everything:
however:
julia> raw.params
Dict{String, Any} with 23 entries:
"protocolName" => "NoName"
"institutionName" => "Pontificia Universidad Catolica de Chile"
"reconSize" => [32, 32, 1]
"enc_lim_repetition" => Limit(0, 0, 0)
"enc_lim_set" => Limit(0, 0, 0)
"enc_lim_segment" => Limit(0, 0, 0)
"userParameters" => Dict{String, Any}("Nblocks"=>2264, "gpu"=>1, "gpu_device"=>0, "type_sim_parts"=>Bool[0, 1, 0, 1, 0, 1, 0, 1, 0, 1 … 1, 0, 1, 0, 1, 0, 1, 0, 1, 0], "precision…
"enc_lim_phase" => Limit(0, 0, 0)
"enc_lim_average" => Limit(0, 0, 0)
"enc_lim_slice" => Limit(0, 0, 0)
"reconFOV" => Float32[258.065, 258.065, 1.0]
"systemModel" => "v0.8.1"
"H1resonanceFrequency_Hz" => 63866203
"patientName" => "brain3D"
"enc_lim_contrast" => Limit(0, 0, 0)
"trajectory" => "other"
"systemFieldStrength_T" => 1.5
"enc_lim_kspace_encoding_step_0" => Limit(0, 31, 16)
"enc_lim_kspace_encoding_step_1" => Limit(0, 31, 16)
"encodedFOV" => Float32[258.065, 258.065, 1.0]
"enc_lim_kspace_encoding_step_2" => Limit(0, 0, 0)
"systemVendor" => "KomaMRI.jl"
"encodedSize" => [32, 32, 1]
with "encodedSize" => [32, 32, 1]
it is not taking the values from seq.DEF as far as I can tell. and the acq parameters get generated correspondingly:
julia> acq = AcquisitionData( raw);
julia> acq.encodingSize
(32, 32)
julia> acq.fov
(258.06451416015625, 258.0649108886719, 1.0)
Trying to manually set the recon size fails with:
julia> Nx, Ny, Nz = [32, 32, 32];
julia> reconParams = Dict{Symbol,Any}(:reco=>"direct", :reconSize=>(Nx, Ny, Nz));
julia> image = reconstruction(acq, reconParams);
ERROR: ArgumentError: Nodes x have dimension 2 != 3
Next step is to fill in the encoded size correctly in raw.
with
"encodedSize" => [32, 32, 1]
it is not taking the values from seq.DEF as far as I can tell. and the acq parameters get generated correspondingly:
Ok, this should work. So it is an unrelated bug (@beorostica can you create an issue?).
julia> image = reconstruction(acq, reconParams); ERROR: ArgumentError: Nodes x have dimension 2 != 3
Last time I tried using 3D k-spaces (nodes) with MRIReco it always gave an error, so we are just storing 2D k-spaces in the raw data (or at least that is the default ndims=2
of signal_to_raw_data
). This could create this dimension mismatch. Maybe we just need to check if :reconSize
has 3 components to save a 3D k-space, otherwise a 2D k-space. This is the main thing we need to explore to fix this issue.
Manually editing raw.params:
julia> raw.params
Dict{String, Any} with 23 entries:
"protocolName" => "NoName"
"institutionName" => "Pontificia Universidad Catolica de Chile"
"reconSize" => [32, 32, 32]
"enc_lim_repetition" => Limit(0, 0, 0)
"enc_lim_set" => Limit(0, 0, 0)
"enc_lim_segment" => Limit(0, 0, 0)
"userParameters" => Dict{String, Any}("Nblocks"=>2264, "gpu"=>1, "gpu_device"=>0, "type_sim_parts"=>Bool[0, 1, 0, 1, 0, 1, 0, 1, 0, 1 … 1, 0, 1, 0, 1, 0, 1, 0, 1, 0], "precision…
"enc_lim_phase" => Limit(0, 0, 0)
"enc_lim_average" => Limit(0, 0, 0)
"enc_lim_slice" => Limit(0, 0, 0)
"reconFOV" => [258.065, 258.065, 258.065]
"systemModel" => "v0.8.1"
"H1resonanceFrequency_Hz" => 63866203
"patientName" => "brain3D"
"enc_lim_contrast" => Limit(0, 0, 0)
"trajectory" => "other"
"systemFieldStrength_T" => 1.5
"enc_lim_kspace_encoding_step_0" => Limit(0, 31, 16)
"enc_lim_kspace_encoding_step_1" => Limit(0, 31, 16)
"encodedFOV" => [258.065, 258.065, 258.065]
"enc_lim_kspace_encoding_step_2" => Limit(0, 31, 16)
"systemVendor" => "KomaMRI.jl"
"encodedSize" => [32, 32, 32]
gives mixed results:
julia> acq = AcquisitionData( raw);
julia> acq.encodingSize
(32, 32)
julia> acq.fov
(258.065, 258.065, 258.065)
julia> acq.traj
1-element Vector{MRIBase.Trajectory{Float32}}:
MRIBase.Trajectory{Float32}("Custom", Float32[-0.20536378 -0.19211452 … 0.19211386 0.20536311; 0.21198764 0.21198764 … 0.21198764 0.21198764], Float32[0.0, 6.2f-5, 0.000124, 0.000186, 0.000248, 0.00031, 0.000372, 0.000434, 0.000496, 0.000558 … 0.001364, 0.001426, 0.001488, 0.00155, 0.001612, 0.001674, 0.001736, 0.001798, 0.00186, 0.001922], 0.0f0, 0.001f0, 32, 32, 1, false, true)
and unfortunately still does not reconstruct:
julia> Nx, Ny, Nz = raw.params["reconSize"][1:3]
3-element Vector{Int64}:
32
32
32
julia> reconParams = Dict{Symbol,Any}(:reco=>"direct", :reconSize=>(Nx, Ny, Nz));
julia> reconParams
Dict{Symbol, Any} with 2 entries:
:reconSize => (32, 32, 32)
:reco => "direct"
julia> image = reconstruction(acq, reconParams)
ERROR: ArgumentError: Nodes x have dimension 2 != 3
Stacktrace:
[1] initParams(k::Matrix{Float64}, N::Tuple{Int64, Int64, Int64}, dims::UnitRange{Int64}; kargs::Base.Pairs{Symbol, Int64, Tuple{Symbol, Symbol}, NamedTuple{(:m, :σ), Tuple{Int64, Int64}}})
@ NFFT ~/.julia/packages/NFFT/RT2hs/src/precomputation.jl:20
[2] initParams
@ ~/.julia/packages/NFFT/RT2hs/src/precomputation.jl:3 [inlined]
[3] NFFT.NFFTPlan(k::Matrix{Float64}, N::Tuple{Int64, Int64, Int64}; dims::UnitRange{Int64}, fftflags::Nothing, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol, Symbol}, NamedTuple{(:m, :σ), Tuple{Int64, Int64}}})
@ NFFT ~/.julia/packages/NFFT/RT2hs/src/implementation.jl:78
[4] NFFTPlan
@ ~/.julia/packages/NFFT/RT2hs/src/implementation.jl:73 [inlined]
[5] macro expansion
@ ~/.julia/packages/NFFT/RT2hs/src/NFFT.jl:49 [inlined]
[6] macro expansion
@ ./timing.jl:393 [inlined]
[7] #plan_nfft#130
@ ~/.julia/packages/NFFT/RT2hs/src/NFFT.jl:48 [inlined]
[8] plan_nfft
@ ~/.julia/packages/NFFT/RT2hs/src/NFFT.jl:46 [inlined]
[9] #plan_nfft#2
@ ~/.julia/packages/AbstractNFFTs/Xd3qS/src/derived.jl:13 [inlined]
[10] plan_nfft
@ ~/.julia/packages/AbstractNFFTs/Xd3qS/src/derived.jl:13 [inlined]
[11] samplingDensity(acqData::AcquisitionData{Float32, 2}, shape::Tuple{Int64, Int64, Int64})
@ MRIBase ~/.julia/packages/MRIBase/oNfYy/src/Datatypes/AcqData.jl:288
[12] setupDirectReco(acqData::AcquisitionData{Float32, 2}, recoParams::Dict{Symbol, Any})
@ MRIReco ~/src/MRIReco.jl/src/Reconstruction/RecoParameters.jl:19
[13] reconstruction(acqData::AcquisitionData{Float32, 2}, recoParams::Dict{Symbol, Any})
@ MRIReco ~/src/MRIReco.jl/src/Reconstruction/Reconstruction.jl:37
[14] top-level scope
@ REPL[80]:1
probably because of acq.traj and other incorrect parameters or dimensions of trajectory and other related data?
The data is getting there but in the wrong dimensionality/indexing:
julia> acq.kdata
1×32×1 Array{Matrix{ComplexF32}, 3}:
[:, :, 1] =
[-5.97573-6.46294im; -21.4331+0.217442im; … ; 21.9411-4.22948im; -11.1961-2.9045im;;] … [22.9858-1.1932im; 5.25461-2.63809im; … ; 19.4627-5.74995im; -4.43083-8.30492im;;]
julia> size( acq.kdata[1, 1, 1])
(1024, 1)
julia> size( acq.kdata[1, 2, 1])
(1024, 1)
julia> size( acq.kdata[1, 32, 1])
(1024, 1)
@cncastillo @beorostica
Line 448 suggests it may work if the DEF items are placed in the seq file. Line 456 may be the problem when not present in 3D sequence? I will try the Nx-z DEF items in the seq file.
if !haskey(seq.DEF,"Nx")
seq.DEF["Nz"] = Nz #Number of unique RF frequencies, in a 3D acquisition this should not work
This is a test where the Nx-z parameters are in the [DEFINITIONS] block of the seq file:
[DEFINITIONS]
AdcRasterTime 6.25e-05
BlockDurationRaster 1e-05
GradientRasterTime 1e-05
RadiofrequencyRasterTime 1e-06
Nx 32
Ny 32
Nz 32
(base) curt@green:~/src/KomaMRI$ julia
_
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.9.4 (2023-11-14)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
julia> using Revise
julia> using NativeFileDialog, MAT, MRIReco
julia> using KomaMRIBase, KomaMRICore, KomaMRIFiles, KomaMRIPlots
julia> sys = Scanner();
julia> seq = read_seq( seq_file)
[ Info: Loading sequence gre3d_cac.seq ...
Sequence[ τ = 42992.22 ms | blocks: 6444 | ADC: 1024 | GR: 7290 | RF: 1074 | DEF: 15 ]
julia> seq.DEF
Dict{String, Any} with 15 entries:
"Nz" => 32.0
"Nx" => 32.0
"signature" => ""
"Ny" => 32.0
"AdcRasterTime" => 1.0e-7
"BlockDurationRaster" => 1.0e-5
"extension" => [0, 0, 0, 0, 0, 0, 0, 0, 0, 0 … 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
"PulseqVersion" => 1004001
"FileName" => "gre3d_cac.seq"
"GradientRasterTime" => 1.0e-5
"TotalDuration" => 42.9922
"FOV" => [0.19, 0.19, 0.19]
"Name" => "gre3d_cac"
"RadiofrequencyRasterTime" => 1.0e-6
"additional_text" => "# Format of extension lists:\n# id type ref next_id\n# next_id of 0 terminates the list\n# Extension list is followed by extension spec…
julia> get_flip_angles( seq[1])
1-element Vector{Float64}:
8.0
julia> raw = simulate(obj, seq, sys, sim_params=sim_params);
┌ Info: Running simulation in the GPU (Quadro RTX 8000)
│ koma_version = v"0.8.1"
│ sim_method = Bloch()
│ spins = 78126
│ time_points = 74945
└ adc_points = 32768
Progress: 100%|██████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:10
simulated_blocks: 2168
rf_blocks: 1075
acq_samples: 32768
10.773805 seconds (9.39 M allocations: 498.809 MiB, 3.71% gc time)
julia> raw.params
Dict{String, Any} with 23 entries:
"protocolName" => "gre3d_cac"
"institutionName" => "Pontificia Universidad Catolica de Chile"
"reconSize" => [32, 32, 1]
"enc_lim_repetition" => Limit(0, 0, 0)
"enc_lim_set" => Limit(0, 0, 0)
"enc_lim_segment" => Limit(0, 0, 0)
"userParameters" => Dict{String, Any}("Nblocks"=>2168, "gpu"=>1, "gpu_device"=>0, "type_sim_parts"=>Bool[0, 1, 0, 1, 0, 1, 0, 1, 0, 1 … 1, 0, 1, 0, …
"enc_lim_phase" => Limit(0, 0, 0)
"enc_lim_average" => Limit(0, 0, 0)
"enc_lim_slice" => Limit(0, 0, 0)
"reconFOV" => Float32[190.0, 190.0, 1.0]
"systemModel" => "v0.8.1"
"H1resonanceFrequency_Hz" => 63866203
"patientName" => "brain3D"
"enc_lim_contrast" => Limit(0, 0, 0)
"trajectory" => "other"
"systemFieldStrength_T" => 1.5
"enc_lim_kspace_encoding_step_0" => Limit(0, 30, 16)
"enc_lim_kspace_encoding_step_1" => Limit(0, 30, 16)
⋮ => ⋮
julia> acq = AcquisitionData( raw);
julia> reconParams = Dict{Symbol,Any}(:reco=>"direct", :reconSize=>(32, 32, 32));
julia> image = reconstruction(acq, reconParams);
ERROR: ArgumentError: Nodes x have dimension 2 != 3
Stacktrace:
[1] initParams(k::Matrix{Float64}, N::Tuple{Int64, Int64, Int64}, dims::UnitRange{Int64}; kargs::Base.Pairs{Symbol, Int64, Tuple{Symbol, Symbol}, NamedTuple{(:m, :σ), Tuple{Int64, Int64}}})
@ NFFT ~/.julia/packages/NFFT/RT2hs/src/precomputation.jl:20
[2] initParams
@ ~/.julia/packages/NFFT/RT2hs/src/precomputation.jl:3 [inlined]
[3] NFFT.NFFTPlan(k::Matrix{Float64}, N::Tuple{Int64, Int64, Int64}; dims::UnitRange{Int64}, fftflags::Nothing, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol, Symbol}, NamedTuple{(:m, :σ), Tuple{Int64, Int64}}})
@ NFFT ~/.julia/packages/NFFT/RT2hs/src/implementation.jl:78
[4] NFFTPlan
@ ~/.julia/packages/NFFT/RT2hs/src/implementation.jl:73 [inlined]
[5] macro expansion
@ ~/.julia/packages/NFFT/RT2hs/src/NFFT.jl:49 [inlined]
[6] macro expansion
@ ./timing.jl:393 [inlined]
[7] #plan_nfft#130
@ ~/.julia/packages/NFFT/RT2hs/src/NFFT.jl:48 [inlined]
[8] plan_nfft
@ ~/.julia/packages/NFFT/RT2hs/src/NFFT.jl:46 [inlined]
[9] #plan_nfft#2
@ ~/.julia/packages/AbstractNFFTs/Xd3qS/src/derived.jl:13 [inlined]
[10] plan_nfft
@ ~/.julia/packages/AbstractNFFTs/Xd3qS/src/derived.jl:13 [inlined]
[11] samplingDensity(acqData::AcquisitionData{Float32, 2}, shape::Tuple{Int64, Int64, Int64})
@ MRIBase ~/.julia/packages/MRIBase/oNfYy/src/Datatypes/AcqData.jl:288
[12] setupDirectReco(acqData::AcquisitionData{Float32, 2}, recoParams::Dict{Symbol, Any})
@ MRIReco ~/src/MRIReco.jl/src/Reconstruction/RecoParameters.jl:19
[13] reconstruction(acqData::AcquisitionData{Float32, 2}, recoParams::Dict{Symbol, Any})
@ MRIReco ~/src/MRIReco.jl/src/Reconstruction/Reconstruction.jl:37
[14] top-level scope
@ REPL[24]:1
Here is some more testing:
julia> raw.params["encodedFOV"]
3-element Vector{Float32}:
190.0
190.0
1.0
julia> raw.params["encodedFOV"]=[190.0, 190.0, 190.0]
3-element Vector{Float64}:
190.0
190.0
190.0
julia> raw.params["encodedSize"]=[32, 32, 32]
3-element Vector{Int64}:
32
32
32
julia> raw.params
Dict{String, Any} with 23 entries:
"protocolName" => "gre3d_cac"
"institutionName" => "Pontificia Universidad Catolica de Chile"
"reconSize" => [32, 32, 32]
"enc_lim_repetition" => Limit(0, 0, 0)
"enc_lim_set" => Limit(0, 0, 0)
"enc_lim_segment" => Limit(0, 0, 0)
"userParameters" => Dict{String, Any}("Nblocks"=>2168, "gpu"=>1, "gpu_device"=>0, "type_sim_parts"=>Bool[0, 1, 0, 1, 0, 1, 0, 1, 0, 1 … 1, 0, 1, 0, 1, 0, 1, 0, 1, 0…
"enc_lim_phase" => Limit(0, 0, 0)
"enc_lim_average" => Limit(0, 0, 0)
"enc_lim_slice" => Limit(0, 0, 0)
"reconFOV" => [190.0, 190.0, 190.0]
"systemModel" => "v0.8.1"
"H1resonanceFrequency_Hz" => 63866203
"patientName" => "brain3D"
"enc_lim_contrast" => Limit(0, 0, 0)
"trajectory" => "other"
"systemFieldStrength_T" => 1.5
"enc_lim_kspace_encoding_step_0" => Limit(0, 30, 16)
"enc_lim_kspace_encoding_step_1" => Limit(0, 30, 16)
"encodedFOV" => [190.0, 190.0, 190.0]
"enc_lim_kspace_encoding_step_2" => Limit(0, 30, 16)
"systemVendor" => "KomaMRI.jl"
"encodedSize" => [32, 32, 32]
julia> acq = AcquisitionData( raw);
julia> reconParams = Dict{Symbol,Any}(:reco=>"direct", :reconSize=>(32, 32, 32));
julia> image = reconstruction(acq, reconParams);
ERROR: ArgumentError: Nodes x have dimension 2 != 3
Stacktrace:
[1] initParams(k::Matrix{Float64}, N::Tuple{Int64, Int64, Int64}, dims::UnitRange{Int64}; kargs::Base.Pairs{Symbol, Int64, Tuple{Symbol, Symbol}, NamedTuple{(:m, :σ), Tuple{Int64, Int64}}})
@ NFFT ~/.julia/packages/NFFT/RT2hs/src/precomputation.jl:20
[2] initParams
@ ~/.julia/packages/NFFT/RT2hs/src/precomputation.jl:3 [inlined]
[3] NFFT.NFFTPlan(k::Matrix{Float64}, N::Tuple{Int64, Int64, Int64}; dims::UnitRange{Int64}, fftflags::Nothing, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol, Symbol}, NamedTuple{(:m, :σ), Tuple{Int64, Int64}}})
@ NFFT ~/.julia/packages/NFFT/RT2hs/src/implementation.jl:78
[4] NFFTPlan
@ ~/.julia/packages/NFFT/RT2hs/src/implementation.jl:73 [inlined]
[5] macro expansion
@ ~/.julia/packages/NFFT/RT2hs/src/NFFT.jl:49 [inlined]
[6] macro expansion
@ ./timing.jl:393 [inlined]
[7] #plan_nfft#130
@ ~/.julia/packages/NFFT/RT2hs/src/NFFT.jl:48 [inlined]
[8] plan_nfft
@ ~/.julia/packages/NFFT/RT2hs/src/NFFT.jl:46 [inlined]
[9] #plan_nfft#2
@ ~/.julia/packages/AbstractNFFTs/Xd3qS/src/derived.jl:13 [inlined]
[10] plan_nfft
@ ~/.julia/packages/AbstractNFFTs/Xd3qS/src/derived.jl:13 [inlined]
[11] samplingDensity(acqData::AcquisitionData{Float32, 2}, shape::Tuple{Int64, Int64, Int64})
@ MRIBase ~/.julia/packages/MRIBase/oNfYy/src/Datatypes/AcqData.jl:288
[12] setupDirectReco(acqData::AcquisitionData{Float32, 2}, recoParams::Dict{Symbol, Any})
@ MRIReco ~/src/MRIReco.jl/src/Reconstruction/RecoParameters.jl:19
[13] reconstruction(acqData::AcquisitionData{Float32, 2}, recoParams::Dict{Symbol, Any})
@ MRIReco ~/src/MRIReco.jl/src/Reconstruction/Reconstruction.jl:37
[14] top-level scope
@ REPL[39]:1
julia> acq.
encodingSize fov kdata sequenceInfo subsampleIndices traj
julia> size(acq.kdata)
(1, 32, 1)
julia> acq.kdata[1,1,1]
1024×1 Matrix{ComplexF32}:
-20.856276f0 + 7.8534117f0im
-23.674004f0 + 0.69720984f0im
-31.270561f0 - 11.370958f0im
-24.04995f0 + 42.965137f0im
39.66914f0 + 48.67775f0im
57.792076f0 + 25.581245f0im
21.598516f0 - 4.7730045f0im
-82.14332f0 - 29.625975f0im
-40.8282f0 - 20.949701f0im
13.959399f0 + 4.1346235f0im
-11.851673f0 + 28.035168f0im
24.585815f0 + 9.329832f0im
56.90302f0 + 7.7491198f0im
17.272808f0 + 13.794415f0im
-1.95158f0 - 26.888247f0im
8.668823f0 - 80.32206f0im
⋮
-20.988113f0 - 13.841843f0im
28.07649f0 - 22.968506f0im
45.602043f0 + 26.476562f0im
55.74817f0 + 22.142883f0im
14.360283f0 + 28.204256f0im
17.00085f0 + 1.0903955f0im
5.857233f0 - 23.596966f0im
-25.678976f0 + 40.5332f0im
20.144512f0 + 9.556927f0im
15.072937f0 - 12.6632595f0im
9.912712f0 + 7.5607147f0im
-23.415226f0 - 46.04351f0im
6.5762672f0 + 27.526182f0im
-19.309336f0 + 33.920307f0im
-28.33278f0 - 93.11535f0im
55.864494f0 + 14.333662f0im
julia> acq.kdata[1,32,1]
1024×1 Matrix{ComplexF32}:
-23.373734f0 + 28.150126f0im
34.316868f0 + 29.423656f0im
1.5463567f0 - 3.6292593f0im
-52.664368f0 - 23.4733f0im
-17.954971f0 - 32.63766f0im
7.2143764f0 - 14.636225f0im
-15.502787f0 + 27.550005f0im
-10.73413f0 + 2.437481f0im
-19.018562f0 + 26.472218f0im
-18.13644f0 + 36.484882f0im
8.976894f0 + 4.1743774f0im
41.983356f0 + 26.89516f0im
46.48584f0 - 12.104265f0im
41.22889f0 - 8.658623f0im
40.798172f0 + 39.47073f0im
5.5201626f0 + 19.247066f0im
⋮
-2.829985f0 - 8.459925f0im
32.78006f0 - 21.172215f0im
38.438377f0 + 26.038582f0im
58.80835f0 + 18.670595f0im
20.287176f0 + 25.217157f0im
14.318125f0 + 8.467522f0im
16.699306f0 - 23.905882f0im
-20.987934f0 + 35.578873f0im
15.35615f0 + 21.750748f0im
12.56637f0 - 13.1215515f0im
14.982704f0 + 2.081315f0im
-32.10189f0 - 36.650608f0im
-4.2777576f0 + 27.28927f0im
-5.1820927f0 + 39.759125f0im
-43.646645f0 - 87.57974f0im
39.00147f0 + 13.069844f0im
julia> size( acq.kdata[1,32,1])
(1024, 1)
It looks like fixing the parameters is not enough, the kdata is still packed up wrong?
Testing 3D non-Cartesian, apparently has the same issues as 3D Cartesian.
(base) curt@green:~/src/KomaMRI$ julia
_
_ _ _(_)_ | Documentation: https://docs.julialang.org
(_) | (_) (_) |
_ _ _| |_ __ _ | Type "?" for help, "]?" for Pkg help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 1.9.4 (2023-11-14)
_/ |\__'_|_|_|\__'_| | Official https://julialang.org/ release
|__/ |
julia> using Markdown
julia> using InteractiveUtils
julia> using Revise
julia> using NativeFileDialog, MAT, MRIReco
julia> using KomaMRIBase, KomaMRICore, KomaMRIFiles, KomaMRIPlots
julia> sys = Scanner()
Scanner
B0: Float64 1.5
B1: Float64 1.0e-5
Gmax: Float64 0.06
Smax: Int64 500
ADC_Δt: Float64 2.0e-6
seq_Δt: Float64 1.0e-5
GR_Δt: Float64 1.0e-5
RF_Δt: Float64 1.0e-6
RF_ring_down_T: Float64 2.0e-5
RF_dead_time_T: Float64 0.0001
ADC_dead_time_T: Float64 1.0e-5
julia> sys.RF_dead_time_T=1e-5
1.0e-5
julia> sys.RF_ring_down_T=10e-6
1.0e-5
julia> sys
Scanner
B0: Float64 1.5
B1: Float64 1.0e-5
Gmax: Float64 0.06
Smax: Int64 500
ADC_Δt: Float64 2.0e-6
seq_Δt: Float64 1.0e-5
GR_Δt: Float64 1.0e-5
RF_Δt: Float64 1.0e-6
RF_ring_down_T: Float64 1.0e-5
RF_dead_time_T: Float64 1.0e-5
ADC_dead_time_T: Float64 1.0e-5
julia> seq_file = pick_file( "/home/curt/src"; filterlist="seq")
Gtk-Message: 19:49:26.593: Failed to load module "canberra-gtk-module"
Gtk-Message: 19:49:26.593: Failed to load module "canberra-gtk-module"
(julia:37505): GLib-GIO-WARNING **: 19:49:26.671: Failed to create file monitor for /home/curt/.config/glib-2.0/settings/keyfile: Unable to find default local file monitor type
"/home/curt/src/pulseq_champaign_imaging_llc/matlab/pulseq_m/zte_petra_cac.seq"
julia> seq = read_seq( seq_file)
[ Info: Loading sequence zte_petra_cac.seq ...
Sequence[ τ = 4328.22 ms | blocks: 1052 | ADC: 515 | GR: 3145 | RF: 536 | DEF: 16 ]
julia> seq.DEF
Dict{String, Any} with 16 entries:
"signature" => "e4359a04de6485703a97f6d7c56107ae"
"AdcRasterTime" => 1.0e-7
"GradientRasterTime" => 1.0e-5
"TotalDuration" => 4.32822
"FOV" => [0.256, 0.256, 0.256]
"Name" => "petra"
"Nz" => 1
"Nx" => 200
"Ny" => 515
"SamplesPerShell" => 515.0
"PulseqVersion" => 1004001
"BlockDurationRaster" => 1.0e-5
"extension" => [0, 0, 0, 0, 0, 0, 0, 0, 0, 0 … 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
"FileName" => "zte_petra_cac.seq"
"additional_text" => "# Format of extension lists:\n# id type ref next_id\n# next_id of 0 terminates the list\n# Extension list is followed by extension specifications\n[EXTENSIONS]\n\n# Exten…
"RadiofrequencyRasterTime" => 1.0e-6
julia> seq_plot = plot_seq(seq; range=[0,100], darkmode=true, slider=true,)
[ Info: Listening on: 127.0.0.1:7953, thread id: 1
julia> get_flip_angles( seq[3])
1-element Vector{Float64}:
3.996
julia> sim_params = KomaMRICore.default_sim_params()
Dict{String, Any} with 9 entries:
"return_type" => "raw"
"Nblocks" => 20
"gpu" => true
"Nthreads" => 1
"gpu_device" => 0
"sim_method" => Bloch()
"precision" => "f32"
"Δt" => 0.001
"Δt_rf" => 5.0e-5
julia> sim_params.
age count idxfloor keys maxprobe ndel slots vals
julia> sim_params["Δt_rf"]
5.0e-5
julia> sim_params["Δt_rf"]=1e-6
1.0e-6
julia> obj = brain_phantom3D(; ss=4, start_end=[178, 182]);
julia> sim_params
Dict{String, Any} with 9 entries:
"return_type" => "raw"
"Nblocks" => 20
"gpu" => true
"Nthreads" => 1
"gpu_device" => 0
"sim_method" => Bloch()
"precision" => "f32"
"Δt" => 0.001
"Δt_rf" => 1.0e-6
julia> raw = simulate(obj, seq, sys, sim_params=sim_params);
┌ Info: Running simulation in the GPU (Quadro RTX 8000)
│ koma_version = v"0.8.1"
│ sim_method = Bloch()
│ spins = 13027
│ time_points = 117639
└ adc_points = 103000
Progress: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| Time: 0:00:13
simulated_blocks: 1091
rf_blocks: 536
acq_samples: 103000
13.879795 seconds (9.95 M allocations: 548.532 MiB, 0.90% gc time)
julia> raw.params
Dict{String, Any} with 23 entries:
"protocolName" => "petra"
"institutionName" => "Pontificia Universidad Catolica de Chile"
"reconSize" => [26, 26, 1]
"enc_lim_repetition" => Limit(0, 0, 0)
"enc_lim_set" => Limit(0, 0, 0)
"enc_lim_segment" => Limit(0, 0, 0)
"userParameters" => Dict{String, Any}("Nblocks"=>1091, "gpu"=>1, "gpu_device"=>0, "type_sim_parts"=>Bool[0, 1, 0, 1, 0, 1, 0, 1, 0, 1 … 1, 0, 1, 0, 1, 0, 1, 0, 1, 0], "precision"=>"f3…
"enc_lim_phase" => Limit(0, 0, 0)
"enc_lim_average" => Limit(0, 0, 0)
"enc_lim_slice" => Limit(0, 0, 0)
"reconFOV" => Float32[256.0, 256.0, 1.0]
"systemModel" => "v0.8.1"
"H1resonanceFrequency_Hz" => 63866203
"patientName" => "brain3D"
"enc_lim_contrast" => Limit(0, 0, 0)
"trajectory" => "other"
"systemFieldStrength_T" => 1.5
"enc_lim_kspace_encoding_step_0" => Limit(0, 24, 13)
"enc_lim_kspace_encoding_step_1" => Limit(0, 25, 13)
"encodedFOV" => Float32[256.0, 256.0, 1.0]
"enc_lim_kspace_encoding_step_2" => Limit(0, 0, 0)
"systemVendor" => "KomaMRI.jl"
"encodedSize" => [25, 26, 1]
julia> acq = AcquisitionData( raw);
julia> begin
Nx, Ny = raw.params["reconSize"][1:2];
reconParams = Dict{Symbol,Any}(:reco=>"direct", :reconSize=>(Nx, Ny));
image = reconstruction(acq, reconParams);;
end
6-dimensional AxisArray{ComplexF32,6,...} with axes:
:x, (0.0:10.24:256.0) mm
:y, (0.0:9.846153846153847:246.15384615384616) mm
:z, (0.0:1.0:0.0) mm
:echos, 1:1
:coils, 1:1
:repetitions, 1:1
And data, a 26×26×1×1×1×1 Array{ComplexF32, 6}:
[:, :, 1, 1, 1, 1] =
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im -0.0236505+0.0992867im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im -0.00777292-0.0254325im 0.012173+0.253888im 0.115967+0.127116im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im -0.0793708+0.124642im 0.239155-0.130121im … 0.0784993+0.198522im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im -0.0569608+0.0734885im -0.132711+0.00239755im -0.00624282+0.0227152im -0.401928-0.0373109im -0.136928+0.164801im 0.0+0.0im 0.0+0.0im
-0.0222453+0.124349im 0.14743+0.162326im 0.0811325+0.282061im -0.248402+0.172537im 0.51513+0.188474im 0.265433+0.066368im 0.257239+0.142123im 0.0+0.0im
-0.12643-0.103366im -0.0969422-0.0783985im 0.0239488+0.0854764im 0.0508624-0.0679528im 0.0785211+0.757174im -0.132184-0.129531im -0.0516888-0.192597im 0.0+0.0im
0.0145846+0.0472643im -0.0715923-0.0530962im 0.0238834+0.221171im -0.228209+0.151941im -0.311803+1.07776im 0.0541115+0.721079im 0.0380341-0.126004im 0.0+0.0im
-0.00197875+0.0977978im -0.0290478+0.0666364im 0.0956486+0.0605195im -0.020544+0.442739im … -0.0646548+1.2661im -0.303796+1.2412im -0.133189-0.060337im 0.0+0.0im
0.037106+0.107172im 0.184824+0.0822097im -0.133251-0.186752im 0.141473+0.988213im -0.092847+1.7226im -0.342548+1.21904im 0.0885918+0.442872im 0.0+0.0im
-0.0738069-0.167435im 0.0514082+0.150021im 0.0903598-0.265952im 0.0595301+0.800169im 0.033736+2.22846im 0.0331009+1.57372im 0.261175+0.729681im 0.0829206-0.137986im
0.129208+0.00561913im -0.171928+0.324073im 0.142662-0.406718im -0.24646+1.12177im 0.0174301+2.41918im -0.0122723+2.06395im 0.196475+0.885627im 0.0+0.0im
-0.194472-0.135396im -0.157046+0.263632im -0.23762-0.0996091im -0.270387+1.47725im -0.029394+2.2901im -0.0107559+1.50242im 0.0865599+0.541747im 0.0+0.0im
0.079245-0.0244986im 0.206436+0.0963664im 0.0358069-0.0455092im 0.218636+0.809432im … 0.168469+2.19265im -0.274536+1.46205im -0.0206982-0.0357301im 0.0+0.0im
0.111329-0.0116827im 0.0941668+0.186362im 0.0414217-0.0693443im 0.261088+0.25694im 0.0234841+1.6421im -0.298419+1.49593im 0.0610674+0.0813188im 0.0+0.0im
-0.099452+0.189673im 0.0295172+0.201996im -0.0141953+0.182957im 0.151031+0.0302758im -0.117989+1.24655im 0.0375125+0.718611im 0.0532909-0.0469774im 0.0+0.0im
0.0+0.0im -0.132238-0.202624im -0.0268976+0.13931im 0.0730343-0.285966im -0.00140481+0.758735im 0.0919741+0.324555im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0678213-0.0374102im 0.077156+0.133285im -0.115163-0.118061im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0567113-0.04648im 0.0242621+0.280348im … -0.061529+0.135817im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0612392+0.0382043im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im … 0.0+0.0im 0.0+0.0im 0.0+0.0im 0.0+0.0im
julia> image_plot = plot_image( abs.( image[:, :, 1]);)
(this does give a 3d to 2d projection image)
julia> reconParams = Dict{Symbol,Any}(:reco=>"direct", :reconSize=>(32, 32, 32));
julia> image = reconstruction(acq, reconParams);
ERROR: ArgumentError: Nodes x have dimension 2 != 3
Stacktrace:
[1] initParams(k::Matrix{Float64}, N::Tuple{Int64, Int64, Int64}, dims::UnitRange{Int64}; kargs::Base.Pairs{Symbol, Int64, Tuple{Symbol, Symbol}, NamedTuple{(:m, :σ), Tuple{Int64, Int64}}})
@ NFFT ~/.julia/packages/NFFT/RT2hs/src/precomputation.jl:20
[2] initParams
@ ~/.julia/packages/NFFT/RT2hs/src/precomputation.jl:3 [inlined]
[3] NFFT.NFFTPlan(k::Matrix{Float64}, N::Tuple{Int64, Int64, Int64}; dims::UnitRange{Int64}, fftflags::Nothing, kwargs::Base.Pairs{Symbol, Int64, Tuple{Symbol, Symbol}, NamedTuple{(:m, :σ), Tuple{Int64, Int64}}})
@ NFFT ~/.julia/packages/NFFT/RT2hs/src/implementation.jl:78
[4] NFFTPlan
@ ~/.julia/packages/NFFT/RT2hs/src/implementation.jl:73 [inlined]
[5] macro expansion
@ ~/.julia/packages/NFFT/RT2hs/src/NFFT.jl:49 [inlined]
[6] macro expansion
@ ./timing.jl:393 [inlined]
[7] #plan_nfft#130
@ ~/.julia/packages/NFFT/RT2hs/src/NFFT.jl:48 [inlined]
[8] plan_nfft
@ ~/.julia/packages/NFFT/RT2hs/src/NFFT.jl:46 [inlined]
[9] #plan_nfft#2
@ ~/.julia/packages/AbstractNFFTs/Xd3qS/src/derived.jl:13 [inlined]
[10] plan_nfft
@ ~/.julia/packages/AbstractNFFTs/Xd3qS/src/derived.jl:13 [inlined]
[11] samplingDensity(acqData::AcquisitionData{Float32, 2}, shape::Tuple{Int64, Int64, Int64})
@ MRIBase ~/.julia/packages/MRIBase/oNfYy/src/Datatypes/AcqData.jl:288
[12] setupDirectReco(acqData::AcquisitionData{Float32, 2}, recoParams::Dict{Symbol, Any})
@ MRIReco ~/src/MRIReco.jl/src/Reconstruction/RecoParameters.jl:19
[13] reconstruction(acqData::AcquisitionData{Float32, 2}, recoParams::Dict{Symbol, Any})
@ MRIReco ~/src/MRIReco.jl/src/Reconstruction/Reconstruction.jl:37
[14] top-level scope
@ REPL[32]:1
@cncastillo @beorostica
Last time I tried using 3D k-spaces (nodes) with MRIReco it always gave an error, so we are just storing 2D k-spaces in the raw data (or at least that is the default
ndims=2
ofsignal_to_raw_data
). This could create this dimension mismatch. Maybe we just need to check if:reconSize
has 3 components to save a 3D k-space, otherwise a 2D k-space. This is the main thing we need to explore to fix this issue.
KomaMRICore/src/rawdata/ISMRMRD.jl definitely is only doing 2D.
function signal_to_raw_data(
signal, seq;
phantom_name="Phantom", sys=Scanner(), sim_params=Dict{String,Any}(), ndims=2
)
version = string(VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "..", "Project.toml"))["version"]))
#Number of samples and FOV
_, ktraj = get_kspace(seq) #kspace information
mink = minimum(ktraj, dims=1)
maxk = maximum(ktraj, dims=1)
Wk = maxk .- mink
Δx = 1 ./ Wk[1:2] #[m] Only x-y
Nx = get(seq.DEF, "Nx", 1)
Ny = get(seq.DEF, "Ny", 1)
Nz = get(seq.DEF, "Nz", 1)
if haskey(seq.DEF, "FOV")
FOVx, FOVy, _ = seq.DEF["FOV"] #[m]
if FOVx > 1 FOVx *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm
if FOVy > 1 FOVy *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm
Nx = round(Int64, FOVx / Δx[1])
Ny = round(Int64, FOVy / Δx[2])
else
FOVx = Nx * Δx[1]
FOVy = Ny * Δx[2]
end
Is the idea to have it do 3D again, or fool it at least?
Functions so far that handle 2D only at the moment.
setup_raw()
handle(w, "recon")
function signal_to_raw_data(
signal, seq;
phantom_name="Phantom", sys=Scanner(), sim_params=Dict{String,Any}(), ndims=2
)
What happened?
When using a 3D pulse sequence (MPRAGE) with brain_phantom3D, the recostruction is not performed propperly. Looks like a 3D k-space is considered as a set of 2D k-spaces or something like this. I attach the .seq file as an example (actually, I couldn't attach it as it is, so you will need to remove "txt" after you save it) MPRAGE_brain_25cmfov.seq.txt
Environment