Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

3D recon #324

Draft
wants to merge 55 commits into
base: master
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
55 commits
Select commit Hold shift + click to select a range
36408a2
3D_recon work in progress.
curtcorum Mar 13, 2024
674c0f2
Line 123 setup_raw() needs modifications for 3D?
curtcorum Mar 13, 2024
d3ff944
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Mar 24, 2024
826fa5e
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Mar 27, 2024
6368e89
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Mar 30, 2024
be10215
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Apr 5, 2024
b1ce21c
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Apr 8, 2024
a7fd74e
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Apr 8, 2024
c6c9e21
Merge branch 'master' into 3D_recon
cncastillo May 15, 2024
8eb57a9
Merge branch 'master' into 3D_recon
curtcorum Jun 4, 2024
f7a7c0b
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Jun 5, 2024
440402d
Doctring updated for Sequence().
curtcorum Jun 6, 2024
971bca4
Fixed typo in Sequence() docstring.
curtcorum Jun 6, 2024
43f01dd
Now exporting default_sim_params.
curtcorum Jun 6, 2024
c669ecc
Changes to signal_to_raw_data, now separate counters for z and slices.
curtcorum Jun 7, 2024
032c3f2
Fixed counter calculation for nz and ns for signal_to_raw_data.
curtcorum Jun 7, 2024
38fd47a
Debugging counters for signal_to_raw_data. No need to reset the scan_…
curtcorum Jun 7, 2024
aa37f37
Counters fixed. Needs some more testing. Ready to try 3D recon with M…
curtcorum Jun 7, 2024
410c403
Fixing 2D vs 3D in progress, 2D currently broken in KomaUI()!!
curtcorum Jun 8, 2024
b5fc898
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Jun 8, 2024
eb3e450
Still debugging...
curtcorum Jun 8, 2024
bbaa408
Merge branch '3D_recon' of github.com:curtcorum/KomaMRI.jl into 3D_recon
curtcorum Jun 8, 2024
c82a735
Still forces 2D recons someplace...???!!!
curtcorum Jun 8, 2024
570bde7
Debugging reconstruction calls from KomaUI().
curtcorum Jun 8, 2024
3bd4209
Debug raw data conversion for ge3d_cac.seq, something funny...!!!
curtcorum Jun 8, 2024
d977f25
A fair amount of work to estimate proper 3D recon parameters from seq…
curtcorum Jun 10, 2024
a4c387c
Update ISMRMRD.jl
curtcorum Jun 12, 2024
dacc56d
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Jun 12, 2024
bd51a01
Update ISMRMRD.jl
curtcorum Jun 12, 2024
9aac620
Update ISMRMRD.jl
curtcorum Jun 12, 2024
05d862a
Update ISMRMRD.jl
curtcorum Jun 12, 2024
8719f0a
Update ISMRMRD.jl
curtcorum Jun 12, 2024
acc981d
Update ISMRMRD.jl
curtcorum Jun 12, 2024
eb9b789
Debugging signal_to_raw_data
curtcorum Jun 12, 2024
9f25e24
Check on Ny, Ns, Nz == 0 for scan counters.
curtcorum Jun 12, 2024
26c6c1b
Defaulting to `use_ndseq=true` for `signal_to_raw_data`.
curtcorum Jun 13, 2024
32b5dfd
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Jun 17, 2024
8b46d5c
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Jun 18, 2024
df5aff8
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Jun 19, 2024
aba7452
Merge branch 'master' into 3D_recon
curtcorum Jun 20, 2024
8f3fc9f
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Jun 24, 2024
6eb7702
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Jul 8, 2024
585304f
Working on central recon parameter estimation in KomaMRICore/src/rawd…
curtcorum Jul 9, 2024
ffea298
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Jul 15, 2024
af870bd
Working on estimate_seq_recon_dimension, working for 2d multislice, n…
curtcorum Jul 16, 2024
c8e76fb
Merge branch 'master' into 3D_recon
curtcorum Jul 16, 2024
12e0d63
Fixing merge.
curtcorum Jul 17, 2024
b101fe1
Whitespace formatting.
curtcorum Jul 17, 2024
c98d92b
Forced type of Nx, Ny, Nz, Ns read from seq.DEF to be Int64.
curtcorum Jul 17, 2024
18539d1
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Jul 17, 2024
0f317cb
LiearAlgebra added to Project. Check for radial sampling.
curtcorum Jul 18, 2024
3d5c2cc
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Aug 5, 2024
af4f6da
Merge branch 'master_merge' into 3D_recon_merge
curtcorum Sep 1, 2024
5bcff4f
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Sep 27, 2024
f393f4a
Merge branch 'JuliaHealth:master' into 3D_recon
curtcorum Oct 3, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 3 additions & 3 deletions KomaMRIBase/src/datatypes/Sequence.jl
Original file line number Diff line number Diff line change
Expand Up @@ -18,9 +18,9 @@ block. This struct serves as an input for the simulation.
- `ADC`: (`::Array{ADC,1}`) ADC block vector
- `DUR`: (`::Vector`, `[s]`) duration block vector
- `DEF`: (`::Dict{String, Any}`) dictionary with relevant information of the sequence.
Possible keys could be [`"AdcRasterTime"`, `"GradientRasterTime"`, `"Name"`, `"Nz"`,
`"Num_Blocks"`, `"Nx"`, `"Ny"`, `"PulseqVersion"`, `"BlockDurationRaster"`,
`"FileName"`, `"RadiofrequencyRasterTime"`]
Possible keys could be [`"GradientRasterTime"`, `"RadiofrequencyRasterTime"`, `"AdcRasterTime"`,
`"BlockDurationRaster"`, `"Name"`, `"FOV"`, `"TE"`, `"TR"`, `"TotalDuration"`, `"Num_Blocks"`,
`"Nx"`, `"Ny"`, `"Nz"`, `"PulseqVersion"`, `"FileName"`]

# Returns
- `seq`: (`::Sequence`) Sequence struct
Expand Down
1 change: 1 addition & 0 deletions KomaMRICore/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ KomaMRIBase = "d0bc0b20-b151-4d03-b2a4-6ca51751cb9c"
ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca"
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

[weakdeps]
AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
Expand Down
6 changes: 4 additions & 2 deletions KomaMRICore/src/KomaMRICore.jl
Original file line number Diff line number Diff line change
Expand Up @@ -23,10 +23,12 @@ include("simulation/SimulatorCore.jl")
include("simulation/Flow.jl")

# ISMRMRD
export signal_to_raw_data
using LinearAlgebra
#export signal_to_raw_data
export signal_to_raw_data, estimate_seq_recon_dimension # *** CAC 240708 for debugging
# Simulator
export Mag
export simulate, simulate_slice_profile
export simulate, simulate_slice_profile, default_sim_params
# Spinors
export Spinor, Rx, Ry, Rz, Q, Un

Expand Down
242 changes: 197 additions & 45 deletions KomaMRICore/src/rawdata/ISMRMRD.jl
Original file line number Diff line number Diff line change
Expand Up @@ -43,7 +43,7 @@ const ISMRMRD_ACQ_USER6 = 1b64 << ( 62 - 1 )
const ISMRMRD_ACQ_USER7 = 1b64 << ( 63 - 1 )
const ISMRMRD_ACQ_USER8 = 1b64 << ( 64 - 1 )
"""
raw = signal_to_raw_data(signal, seq; phantom_name, sys, sim_params)
raw = signal_to_raw_data(signal, seq; phantom_name, sys, sim_params, ndims=2, use_ndseq=false)

Transforms the raw signal into a RawAcquisitionData struct (nearly equivalent to the ISMRMRD
format) used for reconstruction with MRIReco.
Expand All @@ -56,6 +56,8 @@ format) used for reconstruction with MRIReco.
- `phantom_name`: (`::String`, `="Phantom"`) phantom name
- `sys`: (`::Scanner`, `=Scanner()`) Scanner struct
- `sim_params`: (`::Dict{String, Any}`, `=Dict{String,Any}()`) simulation parameter dictionary
- `ndims` : (`::Integer`, `=2`) number of dimensions of the reconstruction
- `use_ndseq` = (`Bool`, `=false`) attempts to estimate dimension of reconstruction from seq file

# Returns
- `raw`: (`::RawAcquisitionData`) RawAcquisitionData struct
Expand All @@ -70,34 +72,31 @@ julia> sim_params = KomaMRICore.default_sim_params(); sim_params["return_type"]

julia> signal = simulate(obj, seq, sys; sim_params)

julia> raw = signal_to_raw_data(signal, seq)
julia> raw = signal_to_raw_data(signal, seq; ndims=3)

julia> plot_signal(raw)
```
"""
function signal_to_raw_data(
signal, seq;
phantom_name="Phantom", sys=Scanner(), sim_params=Dict{String,Any}(), ndims=2
phantom_name="Phantom", sys=Scanner(), sim_params=Dict{String,Any}(), ndims=2, use_ndseq=true
)

#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])
Nd_seq, Nx, Ny, Nz, Ns, FOVx, FOVy, FOVz, Δx, ktraj = estimate_seq_recon_dimension(seq; sim_params)
if use_ndseq
if ndims != Nd_seq; @warn("Using estimated dimension of $Nd_seq for .mrd file."); ndims = Nd_seq; end
else
if ndims != Nd_seq; @warn("Seqfile is $Nd_seq dimensional but .mrd file will be $ndims."); end
end
if ndims == 2
@info "Creating 2D ISMRMRD file ..."
elseif ndims == 3
@info "Creating 3D ISMRMRD file ..."
else
FOVx = Nx * Δx[1]
FOVy = Ny * Δx[2]
@info("Creating a $ndims dimensional ISMRMRD file...")
end

#It needs to be transposed for the raw data
ktraj = maximum(2*abs.(ktraj[:])) == 0 ? transpose(ktraj) : transpose(ktraj)./ maximum(2*abs.(ktraj[:]))

Expand All @@ -124,17 +123,17 @@ function signal_to_raw_data(
# "trajectoryDescription" => Dict{String, Any}("comment"=>""), #You can put wathever you want here: comment, bandwidth, MaxGradient_G_per_cm, MaxSlewRate_G_per_cm_per_s, interleaves, etc
#encoding
# encodedSpace
"encodedSize" => [Nx, Ny, 1], #encodedSpace>matrixSize
"encodedFOV" => Float32.([FOVx, FOVy, 1e-3]*1e3), #encodedSpace>fieldOfView_mm
"encodedSize" => [Nx, Ny, Nz], #encodedSpace>matrixSize
"encodedFOV" => Float32.([FOVx, FOVy, FOVz]*1e3), #encodedSpace>fieldOfView_mm
# reconSpace
"reconSize" => [Nx+Nx%2, Ny+Ny%2, 1], #reconSpace>matrixSize
"reconFOV" => Float32.([FOVx, FOVy, 1e-3]*1e3), #reconSpace>fieldOfView_mm
"reconSize" => [Nx, Ny, Nz], #reconSpace>matrixSize
"reconFOV" => Float32.([FOVx, FOVy, FOVz]*1e3), #reconSpace>fieldOfView_mm
#encodingLimits
"enc_lim_kspace_encoding_step_0" => Limit(0, Nx-1, ceil(Int, Nx / 2)), #min, max, center, e.g. phase encoding line number
"enc_lim_kspace_encoding_step_1" => Limit(0, Ny-1, ceil(Int, Ny / 2)), #min, max, center, e.g. partition encoding number
"enc_lim_kspace_encoding_step_2" => Limit(0, 0, 0), #min, max, center, e.g. partition encoding number
"enc_lim_kspace_encoding_step_0" => Limit(0, Nx-1, floor(Int, Nx / 2)), #min, max, center, e.g. phase encoding line number
"enc_lim_kspace_encoding_step_1" => Limit(0, Ny-1, floor(Int, Ny / 2)), #min, max, center, e.g. partition encoding number
"enc_lim_kspace_encoding_step_2" => Limit(0, Nz-1, floor(Int, Nz / 2)), #min, max, center, e.g. partition encoding number
"enc_lim_average" => Limit(0, 0, 0), #min, max, center, e.g. signal average number
"enc_lim_slice" => Limit(0, 0, 0), #min, max, center, e.g. imaging slice number
"enc_lim_slice" => Limit(0, Ns-1, floor(Int, Ns / 2)), #min, max, center, e.g. imaging slice number
"enc_lim_contrast" => Limit(0, 0, 0), #min, max, center, e.g. echo number in multi-echo
"enc_lim_phase" => Limit(0, 0, 0), #min, max, center, e.g. cardiac phase number
"enc_lim_repetition" => Limit(0, 0, 0), #min, max, center, e.g. dynamic number for dynamic scanning
Expand All @@ -153,11 +152,17 @@ function signal_to_raw_data(
#Then, we define the Profiles
profiles = Profile[]
t_acq = get_adc_sampling_times(seq)
Nadcs = sum(is_ADC_on.(seq))
NadcsPerImage = floor(Int, Nadcs / Nz)
Nro = sum(is_ADC_on.(seq))
#Nyt = max(Ny, 1); Nst = max(Ns, 1); Nzt = max(Nz, 1) #zero is really one...move to estimate_seq_recon_dimension
NroPerPE1 = ceil(Int, Nro / (Ny*Nz*Ns))
NroPerSlice = ceil(Int, Nro / Ns)
NroPerPE2 = ceil(Int, Nro / Nz)
@debug "Nro, NroPerPE1, NroPerSlice, NroPerPE2 = $Nro, $NroPerPE1, $NroPerSlice, $NroPerPE2"
scan_counter = 0
nz = 0
current = 1
ny = 0 #PE1 counter
ns = 0 #Slice counter
nz = 0 #PE2 counter
current = 1 #used for traj and data partition
for s = seq #Iterate over sequence blocks
if is_ADC_on(s)
Nsamples = s.ADC.N[1]
Expand All @@ -167,9 +172,11 @@ function signal_to_raw_data(
if scan_counter == 0
flag += ISMRMRD_ACQ_FIRST_IN_ENCODE_STEP1
flag += ISMRMRD_ACQ_FIRST_IN_SLICE
elseif scan_counter == Nadcs - 1
if Nz > 1; flag += ISMRMRD_ACQ_FIRST_IN_ENCODE_STEP2; end
elseif scan_counter == Nro - 1 #Needs fixing? CAC 240609 ***
flag += ISMRMRD_ACQ_LAST_IN_ENCODE_STEP1
flag += ISMRMRD_ACQ_LAST_IN_SLICE
if Nz > 1; flag += ISMRMRD_ACQ_LAST_IN_ENCODE_STEP2; end
end
#Header of profile data, head::AcquisitionHeader
head = AcquisitionHeader(
Expand All @@ -195,15 +202,15 @@ function signal_to_raw_data(
Float32.((0, 0, 1)), #slice_dir float32x3: Directional cosines of the slice direction
Float32.((0, 0, 0)), #patient_table_position float32x3: Patient table off-center
EncodingCounters( #idx uint16x17: Encoding loop counters
UInt16(scan_counter), #kspace_encode_step_1 uint16: e.g. phase encoding line number
UInt16(0), #kspace_encode_step_2 uint16: e.g. partition encoding number
UInt16(0), #average uint16: e.g. signal average number
UInt16(nz), #slice uint16: e.g. imaging slice number
UInt16(0), #contrast uint16: e.g. echo number in multi-echo
UInt16(0), #phase uint16: e.g. cardiac phase number
UInt16(0), #repetition uint16: e.g. dynamic number for dynamic scanning
UInt16(0), #set uint16: e.g. flow encoding set
UInt16(0), #segment uint16: e.g. segment number for segmented acquisition
UInt16(ny), #kspace_encode_step_1 uint16: e.g. phase encoding line number
UInt16(nz), #kspace_encode_step_2 uint16: e.g. partition encoding number
UInt16(0), #average uint16: e.g. signal average number
UInt16(ns), #slice uint16: e.g. imaging slice number
UInt16(0), #contrast uint16: e.g. echo number in multi-echo
UInt16(0), #phase uint16: e.g. cardiac phase number
UInt16(0), #repetition uint16: e.g. dynamic number for dynamic scanning
UInt16(0), #set uint16: e.g. flow encoding set
UInt16(0), #segment uint16: e.g. segment number for segmented acquisition
Tuple(UInt16(0) for i=1:8) #user uint16x8: Free user parameters
),
Tuple(Int32(0) for i=1:8), #user_int int32x8: Free user parameters
Expand All @@ -216,18 +223,163 @@ function signal_to_raw_data(
#Saving profile
push!(profiles, Profile(head, Float32.(traj), ComplexF32.(dat)))
#Update counters
scan_counter += 1
scan_counter += 1 #another ro
current += Nsamples
if scan_counter % NadcsPerImage == 0 #For now only Nz is considered
if scan_counter % NroPerPE1 == 0
ny += 1 #another kspace_encode_step_1
end
ny = ny % Ny
if scan_counter % NroPerPE2 == 0
nz += 1 #another image
scan_counter = 0 #reset counter
end
nz = nz % Nz
if scan_counter % NroPerSlice == 0
ns += 1 #another slice
end
ns = ns % Ns
# more here for other counters, i.e. dynamic
end
end

@debug "scan_counter, ny, nz, ns = $scan_counter, $ny, $nz, $ns"
return RawAcquisitionData(params, profiles)
end

"""
Nd_seq, Nx, Ny, Nz, Ns, FOVx, FOVy, FOVz, Δx, ktraj = estimate_seq_recon_dimension(seq; sim_params)

Utility function for the best estimate of the reconstruction dimension.

# Arguments
- `seq`: (`::Sequence`) Sequence struct

# Keywords
- `sim_params`: (`::Dict{String, Any}`, `=Dict{String,Any}()`) simulation parameter dictionary (IGNORED)

# Returns
- `Nd_seq`: (`Int`) Estimated reconstruction dimension of seq.

# Examples
```julia-repl
julia> seq_file = joinpath(dirname(pathof(KomaMRI)), "../examples/1.sequences/epi_se.seq")

julia> sys, obj, seq = Scanner(), brain_phantom2D(), read_seq(seq_file)

julia> sim_params = KomaMRICore.default_sim_params(); sim_params["return_type"] = "mat"

julia> Nd_seq = estimate_seq_recon_dimension(seq; sim_params)
```
"""
function estimate_seq_recon_dimension(seq; sim_params=Dict{String,Any}(),
)
#remove signal and sim_params...not needed

#K-space info
_, ktraj = get_kspace(seq) #kspace information
mink = minimum(ktraj, dims=1)
maxk = maximum(ktraj, dims=1)
Wk = maxk .- mink
idxs_zero = findall(iszero, Wk) #check for zeros
Wk_el_iszero = iszero.( Wk) # bool array for later
Wk[idxs_zero] .= 1.0e-6
Δx = 1 ./ Wk[1:3] #[m] x-y-z
k_radius= mapslices(norm, ktraj, dims=2) #limit to first N readouts?
k_radius_norm = k_radius ./ maximum(k_radius)

#estimate sequence acquisition structure
Ns_seq = length(unique(seq.RF.Δf)) #slices, slabs, or Cartesean off-center, preps need to & with ADC *** CAC 240708
Np_seq = maximum(adc.N for adc = seq.ADC) #readout length
Nv_seq = sum(map(is_ADC_on, seq)) #total number of readouts or views
@debug "Ns_seq, Np_seq, Nv_seq = $Ns_seq, $Np_seq, $Nv_seq"

#estimate FOV and N from sequence and K-space info
FOVx_k = Np_seq*Δx[1]
FOVy_k = Np_seq*Δx[2]
FOVz_k = Np_seq*Δx[3]
Nx_k = ceil(Int64, FOVx_k / Δx[1])
Ny_k = ceil(Int64, FOVy_k / Δx[2])
Nz_k = ceil(Int64, FOVz_k / Δx[3])
@debug "FOVx_k, FOVy_k, FOVz_k = $FOVx_k, $FOVy_k, $FOVz_k"
@debug "Nx_k, Ny_k, Nz_k = $Nx_k, $Ny_k, $Nz_k"

# some guesses
seq_2d=false; seq_3d=false; seq_cartesean=false; seq_radial=false; seq_spiral=false
if FOVz_k > 100
seq_2d=true
elseif FOVz_k > 0
seq_3d=true
else
@warn "This seq file does not seem to be an imaging sequence."
end
if Np_seq > 5*Nv_seq
seq_spiral = true
elseif sum( k_radius_norm .< 3/Np_seq) >= Nv_seq #Center of K-space highly sampled
seq_radial=true
else
seq_cartesean=true
end
@debug "seq_2d, seq_3d, seq_cartesean, seq_radial, seq_spiral = $seq_2d, $seq_3d, $seq_cartesean, $seq_radial, $seq_spiral"

# Guessing Cartesean recon dimensions
# ideally all estimates of recon dimensions in one place, as late as possible, non-Cartesean ?? *** CAC 240708
Ns = Int64(get(seq.DEF, "Ns", Ns_seq)) #slices or slabs
Nx = Int64(get(seq.DEF, "Nx", Np_seq))
if seq_cartesean
if seq_2d
Ny = Int64(get(seq.DEF, "Ny", ceil(Int64, Nv_seq/Ns))) #pe1
Nz = 1
elseif seq_3d
Ny = Int64(get(seq.DEF, "Ny", Ny_k)) #pe1
Nz = Int64(get(seq.DEF, "Nz", ceil(Int64, Nv_seq/Ny_k))) #pe2
end
else
@warn "Non-Cartesean recon is still under development."
Ny = ceil(Int64, Nv_seq)
if Ny == 1 Ns = 1 end #sinle shot spiral
Nz = 1
end

# explicit reads from seq.DEF
#Nx = ceil(Int64, get(seq.DEF, "Nx", 1))
#Ny = ceil(Int64, get(seq.DEF, "Ny", 1))
#Nz = ceil(Int64, get(seq.DEF, "Nz", 1))
#Ns = ceil(Int64, get(seq.DEF, "Ns", 1)) # number of slices or slabs

#if Nx == 1 Nx = ceil(Int64, FOVx / Δx[1]) end
#if Ny == 1 Ny = ceil(Int64, FOVy / Δx[2]) end
#if Nz == 1 Nz = ceil(Int64, FOVz / Δx[3]) end

if haskey(seq.DEF, "FOV") #needs info or warning for other substitutions???
FOVx, FOVy, FOVz = seq.DEF["FOV"] #[m]
if FOVx > 1.0 #do this in steps until < 1?? seems like mm or cm possible
FOVx *= 1e-2
FOVy *= 1e-2
FOVz *= 1e-2
@warn "Scaling FOV to m from older pulseq file cm."
end
if seq_3d
FOVz = FOVz
else
FOVz = 0.0
end
elseif seq_cartesean
@warn "Estimating FOV parameters from k-space."
FOVx = FOVx_k
FOVy = FOVy_k
if seq_3d
FOVz = FOVz_k
else
FOVz = 0.0
end
end
Nd_seq = (FOVx > .05) + (FOVy > .05) + (FOVz > .05)
s_ktraj = size( ktraj)
@debug "Nd_seq = $Nd_seq"
@debug "Nx, Ny, Nz, Ns = $Nx, $Ny, $Nz, $Ns"
@debug "FOVx, FOVy, FOVz = $FOVx, $FOVy, $FOVz"
@debug "Δx, size( ktraj) = $Δx, $s_ktraj"
return Nd_seq, Nx, Ny, Nz, Ns, FOVx, FOVy, FOVz, Δx, ktraj
end

Base.show(io::IO, raw::RawAcquisitionData) = begin
Nt, Nc = size(raw.profiles[1].data)
compact = get(io, :compact, false)
Expand Down
3 changes: 1 addition & 2 deletions KomaMRICore/src/simulation/GPUFunctions.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,6 @@ _print_devices(backend) = @error "_print_devices called with invalid backend typ
_print_devices(::KA.CPU) = @info "CPU: $(length(Sys.cpu_info())) x $(Sys.cpu_info()[1].model)"
name(::KA.CPU) = "CPU"
set_device!(backend, val) = @error "set_device! called with invalid parameter types: '$(typeof(backend))', '$(typeof(val))'"
set_device!(val) = set_device!(get_backend(true), val)

#oneAPI.jl doesn't support cis (https://github.com/JuliaGPU/oneAPI.jl/pull/443), so
#for now we use a custom function for each backend to implement
Expand Down Expand Up @@ -97,4 +96,4 @@ function print_devices(use_gpu = true)
_print_devices(backend)
end

export print_devices, set_device!
export print_devices
8 changes: 5 additions & 3 deletions KomaMRIFiles/src/Sequence/Pulseq.jl
Original file line number Diff line number Diff line change
Expand Up @@ -468,9 +468,11 @@ function read_seq(filename)
seq.DEF["PulseqVersion"] = version_combined
seq.DEF["signature"] = signature
# Guessing recon dimensions
seq.DEF["Nx"] = get(seq.DEF, "Nx", maximum(adc.N for adc = seq.ADC))
seq.DEF["Nz"] = get(seq.DEF, "Nz", length(unique(seq.RF.Δf)))
seq.DEF["Ny"] = get(seq.DEF, "Ny", sum(map(is_ADC_on, seq)) ÷ seq.DEF["Nz"])
# ideally all estimates of recon dimensions in one place, as late as possible, i.e. not here, non-Cartesean ?? *** CAC 240708
#seq.DEF["Ns"] = get(seq.DEF, "Ns", length(unique(seq.RF.Δf))) #slices or slabs
#seq.DEF["Nx"] = get(seq.DEF, "Nx", maximum(adc.N for adc = seq.ADC)) #readout length
#seq.DEF["Ny"] = get(seq.DEF, "Ny", sum(map(is_ADC_on, seq)) ÷ seq.DEF["Nx"]) #pe1
#seq.DEF["Nz"] = get(seq.DEF, "Nz", sum(map(is_ADC_on, seq)) ÷ seq.DEF["Ny"]) #pe2
#Koma sequence
return seq
end
Expand Down
2 changes: 2 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -8,9 +8,11 @@ AssetRegistry = "bf4720bc-e11a-5d0c-854e-bdca1663c893"
Blink = "ad839575-38b3-5650-b840-f874b8c74a25"
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
Interact = "c601a237-2ae4-5e1e-952c-7a85b0c7eef1"
KomaMRIBase = "d0bc0b20-b151-4d03-b2a4-6ca51751cb9c"
KomaMRICore = "4baa4f4d-2ae9-40db-8331-a7d1080e3f4e"
KomaMRIFiles = "fcf631a6-1c7e-4e88-9e64-b8888386d9dc"
KomaMRIPlots = "76db0263-63f3-4d26-bb9a-5dba378db904"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MAT = "23992714-dd62-5051-b70f-ba57cb901cac"
MRIReco = "bdf86e05-2d2b-5731-a332-f3fe1f9e047f"
Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
Expand Down
Loading