From 36408a22c767a2f58333d4d204d7b184358766b1 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Tue, 12 Mar 2024 22:21:20 -0500 Subject: [PATCH 01/30] 3D_recon work in progress. modified: ../KomaMRICore/src/rawdata/ISMRMRD.jl modified: KomaUI.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 19 ++++++++++++------- src/KomaUI.jl | 4 ++-- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 4951f5477..7887d4ae0 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -85,19 +85,24 @@ function signal_to_raw_data( mink = minimum(ktraj, dims=1) maxk = maximum(ktraj, dims=1) Wk = maxk .- mink - Δx = 1 ./ Wk[1:2] #[m] Only x-y + Δx = 1 ./ Wk[1:3] #[m] Only x-y Nx = get(seq.DEF, "Nx", 1) Ny = get(seq.DEF, "Ny", 1) Nz = get(seq.DEF, "Nz", 1) + Nd_seq = (Nx > 1) + (Ny > 1) + (Nz > 1) + if ndims != Nd_seq; @warn("Seqfile is $Nd_seq dimensional but recon is $ndims."); end if haskey(seq.DEF, "FOV") - FOVx, FOVy, _ = seq.DEF["FOV"] #[m] + FOVx, FOVy, FOVz = 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 + if FOVz > 1 FOVz *= 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]) + Nz = round(Int64, FOVy / Δx[3]) else FOVx = Nx * Δx[1] FOVy = Ny * Δx[2] + FOVz = Nz * Δx[3] end #It needs to be transposed for the raw data ktraj = maximum(2*abs.(ktraj[:])) == 0 ? transpose(ktraj) : transpose(ktraj)./ maximum(2*abs.(ktraj[:])) @@ -125,15 +130,15 @@ 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+Nx%2, Ny+Ny%2, Nz+Nz%2], #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_2" => Limit(0, Nz-1, ceil(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_contrast" => Limit(0, 0, 0), #min, max, center, e.g. echo number in multi-echo diff --git a/src/KomaUI.jl b/src/KomaUI.jl index f4bd93759..fb52033cd 100644 --- a/src/KomaUI.jl +++ b/src/KomaUI.jl @@ -221,8 +221,8 @@ function KomaUI(; darkmode=true, frame=true, phantom_mode="2D", sim=Dict{String, acq_data = AcquisitionData(raw_aux) acq_data.traj[1].circular = false #Removing circular window acq_data.traj[1].nodes = acq_data.traj[1].nodes[1:2,:] ./ maximum(2*abs.(acq_data.traj[1].nodes[:])) #Normalize k-space to -.5 to .5 for NUFFT - Nx, Ny = raw_aux.params["reconSize"][1:2] - rec_params[:reconSize] = (Nx, Ny) + Nx, Ny, Nz = raw_aux.params["reconSize"][1:3] + rec_params[:reconSize] = (Nx, Ny, Nz) rec_params[:densityWeighting] = true # Perform reconstruction From 674c0f2736f7d633ab4111a8d538f1ca7bbea83d Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Tue, 12 Mar 2024 22:23:50 -0500 Subject: [PATCH 02/30] Line 123 setup_raw() needs modifications for 3D? modified: ExportUIFunctions.jl --- src/ui/ExportUIFunctions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/ui/ExportUIFunctions.jl b/src/ui/ExportUIFunctions.jl index 3afa48f8b..cedd7cb9f 100644 --- a/src/ui/ExportUIFunctions.jl +++ b/src/ui/ExportUIFunctions.jl @@ -120,7 +120,7 @@ end """ Returns the default raw signal used by the UI. """ -function setup_raw() +function setup_raw() # needs modifications for 3D? # Define the default RawAcquisitionData struct raw = RawAcquisitionData( From 440402da0ebd3a00e7aa12fe46765fd0d7b84725 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Thu, 6 Jun 2024 17:09:28 -0500 Subject: [PATCH 03/30] Doctring updated for Sequence(). modified: KomaMRIBase/src/datatypes/Sequence.jl --- KomaMRIBase/src/datatypes/Sequence.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/KomaMRIBase/src/datatypes/Sequence.jl b/KomaMRIBase/src/datatypes/Sequence.jl index 5501040a2..c1791088d 100644 --- a/KomaMRIBase/src/datatypes/Sequence.jl +++ b/KomaMRIBase/src/datatypes/Sequence.jl @@ -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 From 971bca42218a416ced838b4e87f330d426722c1d Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Thu, 6 Jun 2024 17:14:45 -0500 Subject: [PATCH 04/30] Fixed typo in Sequence() docstring. modified: KomaMRIBase/src/datatypes/Sequence.jl --- KomaMRIBase/src/datatypes/Sequence.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/KomaMRIBase/src/datatypes/Sequence.jl b/KomaMRIBase/src/datatypes/Sequence.jl index c1791088d..27017f530 100644 --- a/KomaMRIBase/src/datatypes/Sequence.jl +++ b/KomaMRIBase/src/datatypes/Sequence.jl @@ -19,7 +19,7 @@ block. This struct serves as an input for the simulation. - `DUR`: (`::Vector`, `[s]`) duration block vector - `DEF`: (`::Dict{String, Any}`) dictionary with relevant information of the sequence. Possible keys could be [`"GradientRasterTime"`, `"RadiofrequencyRasterTime"`, `"AdcRasterTime"`, - `"BlockDurationRaster"`, `"Name"`, `"FOV"`, `"TE"`, `"TR'", `"TotalDuration"`, `"Num_Blocks"`, + `"BlockDurationRaster"`, `"Name"`, `"FOV"`, `"TE"`, `"TR"`, `"TotalDuration"`, `"Num_Blocks"`, `"Nx"`, `"Ny"`, `"Nz"`, `"PulseqVersion"`, `"FileName"`] # Returns From 43f01dd6e60a328dd33fa7872506651a0b4ed7f5 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Thu, 6 Jun 2024 18:13:15 -0500 Subject: [PATCH 05/30] Now exporting default_sim_params. modified: KomaMRICore/src/KomaMRICore.jl --- KomaMRICore/src/KomaMRICore.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/KomaMRICore/src/KomaMRICore.jl b/KomaMRICore/src/KomaMRICore.jl index b61f4b125..30e7602a3 100644 --- a/KomaMRICore/src/KomaMRICore.jl +++ b/KomaMRICore/src/KomaMRICore.jl @@ -25,7 +25,7 @@ include("simulation/SimulatorCore.jl") export signal_to_raw_data # Simulator export Mag -export simulate, simulate_slice_profile +export simulate, simulate_slice_profile, default_sim_params # Spinors export Spinor, Rx, Ry, Rz, Q, Un From c669ecce5373782e494e63aed7602d3736fa4176 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Thu, 6 Jun 2024 21:27:59 -0500 Subject: [PATCH 06/30] Changes to signal_to_raw_data, now separate counters for z and slices. modified: KomaMRICore/src/rawdata/ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 23 +++++++++++++++-------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 7887d4ae0..55e8a00f4 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -89,6 +89,7 @@ function signal_to_raw_data( Nx = get(seq.DEF, "Nx", 1) Ny = get(seq.DEF, "Ny", 1) Nz = get(seq.DEF, "Nz", 1) + Ns = get(seq.DEF, "Ns", 1) # number of slices or slabs Nd_seq = (Nx > 1) + (Ny > 1) + (Nz > 1) if ndims != Nd_seq; @warn("Seqfile is $Nd_seq dimensional but recon is $ndims."); end if haskey(seq.DEF, "FOV") @@ -130,17 +131,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, Nz], #encodedSpace>matrixSize + "encodedSize" => [Nx, Ny, Nz], #encodedSpace>matrixSize "encodedFOV" => Float32.([FOVx, FOVy, FOVz]*1e3), #encodedSpace>fieldOfView_mm # reconSpace - "reconSize" => [Nx+Nx%2, Ny+Ny%2, Nz+Nz%2], #reconSpace>matrixSize + "reconSize" => [Nx+Nx%2, Ny+Ny%2, Nz+Nz%2], #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, Nz-1, ceil(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, ceil(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 @@ -160,8 +161,10 @@ function signal_to_raw_data( profiles = Profile[] t_acq = get_adc_sampling_times(seq) Nadcs = sum(is_ADC_on.(seq)) - NadcsPerImage = floor(Int, Nadcs / Nz) + NadcsPerSlice = floor(Int, Nadcs / Ns) + NadcsPerPE1 = floor(Int, Nadcs / Nz) scan_counter = 0 + ns = 0 nz = 0 current = 1 for s = seq #Iterate over sequence blocks @@ -202,9 +205,9 @@ function signal_to_raw_data( 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(nz), #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(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 @@ -224,8 +227,12 @@ function signal_to_raw_data( #Update counters scan_counter += 1 current += Nsamples - if scan_counter % NadcsPerImage == 0 #For now only Nz is considered - nz += 1 #another image + if scan_counter % NadcsPerSlice == 0 + ns += 1 #another slice + scan_counter = 0 #reset counter + end + if scan_counter % NadcsPerPE1 == 0 #Using Ns for slice number + nz += 1 #another kspace_encode_step_2 scan_counter = 0 #reset counter end end From 032c3f2b38faa2d4a81e5852b214a0d788f5837f Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Thu, 6 Jun 2024 21:45:25 -0500 Subject: [PATCH 07/30] Fixed counter calculation for nz and ns for signal_to_raw_data. modified: KomaMRICore/src/rawdata/ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 55e8a00f4..1746ec5b8 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -230,8 +230,7 @@ function signal_to_raw_data( if scan_counter % NadcsPerSlice == 0 ns += 1 #another slice scan_counter = 0 #reset counter - end - if scan_counter % NadcsPerPE1 == 0 #Using Ns for slice number + elseif scan_counter % NadcsPerPE1 == 0 #Using Ns for slice number nz += 1 #another kspace_encode_step_2 scan_counter = 0 #reset counter end From 38fd47a31ee0a04618ca2e563e9b86676debec93 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Thu, 6 Jun 2024 21:58:37 -0500 Subject: [PATCH 08/30] Debugging counters for signal_to_raw_data. No need to reset the scan_counter? modified: KomaMRICore/src/rawdata/ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 1746ec5b8..25bae07ff 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -229,11 +229,13 @@ function signal_to_raw_data( current += Nsamples if scan_counter % NadcsPerSlice == 0 ns += 1 #another slice - scan_counter = 0 #reset counter - elseif scan_counter % NadcsPerPE1 == 0 #Using Ns for slice number + #scan_counter = 0 #reset counter + end + if scan_counter % NadcsPerPE1 == 0 #Using Ns for slice number nz += 1 #another kspace_encode_step_2 - scan_counter = 0 #reset counter + #scan_counter = 0 #reset counter end + # more here for other counters, i.e. dynamic end end From aa37f3778983743bd5049e1fd40a20b28a708987 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Thu, 6 Jun 2024 22:46:11 -0500 Subject: [PATCH 09/30] Counters fixed. Needs some more testing. Ready to try 3D recon with MRIReco. Need to configure reco_parameters. modified: KomaMRICore/src/rawdata/ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 25bae07ff..811532508 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -161,8 +161,8 @@ function signal_to_raw_data( profiles = Profile[] t_acq = get_adc_sampling_times(seq) Nadcs = sum(is_ADC_on.(seq)) - NadcsPerSlice = floor(Int, Nadcs / Ns) - NadcsPerPE1 = floor(Int, Nadcs / Nz) + NadcsPerSlice = floor(Int, Nadcs / Ns) +1 + NadcsPerPE1 = floor(Int, Nadcs / Nz) +1 scan_counter = 0 ns = 0 nz = 0 From 410c403528db7e68919c63ebf48c511a01b77d67 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Fri, 7 Jun 2024 20:46:37 -0500 Subject: [PATCH 10/30] Fixing 2D vs 3D in progress, 2D currently broken in KomaUI()!! modified: ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 811532508..50a46f05b 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -85,7 +85,10 @@ function signal_to_raw_data( mink = minimum(ktraj, dims=1) maxk = maximum(ktraj, dims=1) Wk = maxk .- mink - Δx = 1 ./ Wk[1:3] #[m] Only x-y + idxs_zero = findall(iszero, Wk) #check for zeros + @debug Wk + Wk[idxs_zero] .= 1.0 #replace zero elements + Δx = 1 ./ Wk[1:3] #[m] x-y-z Nx = get(seq.DEF, "Nx", 1) Ny = get(seq.DEF, "Ny", 1) Nz = get(seq.DEF, "Nz", 1) From eb3e4501619df7e97f1aa16feee1b7e1532a7780 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Fri, 7 Jun 2024 23:24:41 -0500 Subject: [PATCH 11/30] Still debugging... modified: ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 50a46f05b..708470b3a 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -86,9 +86,10 @@ function signal_to_raw_data( maxk = maximum(ktraj, dims=1) Wk = maxk .- mink idxs_zero = findall(iszero, Wk) #check for zeros - @debug Wk - Wk[idxs_zero] .= 1.0 #replace zero elements + @debug Wk idxs_zero + Wk[idxs_zero] .= 1e6 Δx = 1 ./ Wk[1:3] #[m] x-y-z + @debug "After fixes:" Wk Δx Nx = get(seq.DEF, "Nx", 1) Ny = get(seq.DEF, "Ny", 1) Nz = get(seq.DEF, "Nz", 1) @@ -102,11 +103,11 @@ function signal_to_raw_data( if FOVz > 1 FOVz *= 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]) - Nz = round(Int64, FOVy / Δx[3]) + Nz = round(Int64, FOVy / Δx[3]) else FOVx = Nx * Δx[1] FOVy = Ny * Δx[2] - FOVz = Nz * Δx[3] + FOVz = Ny * Δx[3] end #It needs to be transposed for the raw data ktraj = maximum(2*abs.(ktraj[:])) == 0 ? transpose(ktraj) : transpose(ktraj)./ maximum(2*abs.(ktraj[:])) @@ -140,11 +141,11 @@ function signal_to_raw_data( "reconSize" => [Nx+Nx%2, Ny+Ny%2, Nz+Nz%2], #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, Nz-1, ceil(Int, Nz / 2)), #min, max, center, e.g. partition encoding number + "enc_lim_kspace_encoding_step_0" => Limit(0, Nx-1, ceil(Int, Nx / 2)-1), #min, max, center, e.g. phase encoding line number + "enc_lim_kspace_encoding_step_1" => Limit(0, Ny-1, ceil(Int, Ny / 2)-1), #min, max, center, e.g. partition encoding number + "enc_lim_kspace_encoding_step_2" => Limit(0, Nz-1, ceil(Int, Nz / 2)-1), #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, Ns-1, ceil(Int, Ns / 2)), #min, max, center, e.g. imaging slice number + "enc_lim_slice" => Limit(0, Ns-1, ceil(Int, Ns / 2)-1), #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 From c82a735142f544a556e76406125af5eac1ed3269 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Sat, 8 Jun 2024 01:35:00 -0500 Subject: [PATCH 12/30] Still forces 2D recons someplace...???!!! modified: KomaMRICore/src/rawdata/ISMRMRD.jl modified: src/KomaUI.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 23 +++++++++++++---------- src/KomaUI.jl | 15 ++++++++++----- 2 files changed, 23 insertions(+), 15 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 708470b3a..bf7b2983b 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -77,7 +77,7 @@ 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=-1, ) version = string(VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "..", "Project.toml"))["version"])) #Number of samples and FOV @@ -86,16 +86,12 @@ function signal_to_raw_data( maxk = maximum(ktraj, dims=1) Wk = maxk .- mink idxs_zero = findall(iszero, Wk) #check for zeros - @debug Wk idxs_zero Wk[idxs_zero] .= 1e6 Δx = 1 ./ Wk[1:3] #[m] x-y-z - @debug "After fixes:" Wk Δx Nx = get(seq.DEF, "Nx", 1) Ny = get(seq.DEF, "Ny", 1) Nz = get(seq.DEF, "Nz", 1) Ns = get(seq.DEF, "Ns", 1) # number of slices or slabs - Nd_seq = (Nx > 1) + (Ny > 1) + (Nz > 1) - if ndims != Nd_seq; @warn("Seqfile is $Nd_seq dimensional but recon is $ndims."); end if haskey(seq.DEF, "FOV") FOVx, FOVy, FOVz = seq.DEF["FOV"] #[m] if FOVx > 1 FOVx *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm @@ -109,6 +105,13 @@ function signal_to_raw_data( FOVy = Ny * Δx[2] FOVz = Ny * Δx[3] end + Nd_seq = (Nx > 1) + (Ny > 1) + (Nz > 1) + if ndims > -1 #just in case ndims is provided in a call + if ndims != Nd_seq; @warn("Seqfile is $Nd_seq dimensional but .mrd file will be $ndims."); end + else + ndims = Nd_seq + end + @info "Creating $ndims dimensional ISMRMRD file ..." #It needs to be transposed for the raw data ktraj = maximum(2*abs.(ktraj[:])) == 0 ? transpose(ktraj) : transpose(ktraj)./ maximum(2*abs.(ktraj[:])) @@ -138,14 +141,14 @@ function signal_to_raw_data( "encodedSize" => [Nx, Ny, Nz], #encodedSpace>matrixSize "encodedFOV" => Float32.([FOVx, FOVy, FOVz]*1e3), #encodedSpace>fieldOfView_mm # reconSpace - "reconSize" => [Nx+Nx%2, Ny+Ny%2, Nz+Nz%2], #reconSpace>matrixSize + "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)-1), #min, max, center, e.g. phase encoding line number - "enc_lim_kspace_encoding_step_1" => Limit(0, Ny-1, ceil(Int, Ny / 2)-1), #min, max, center, e.g. partition encoding number - "enc_lim_kspace_encoding_step_2" => Limit(0, Nz-1, ceil(Int, Nz / 2)-1), #min, max, center, e.g. partition encoding number + "enc_lim_kspace_encoding_step_0" => Limit(0, Nx-1, ceil(Int, Nx / 2)-1),#min, max, center, e.g. phase encoding line number + "enc_lim_kspace_encoding_step_1" => Limit(0, Ny-1, ceil(Int, Ny / 2)-1),#min, max, center, e.g. partition encoding number + "enc_lim_kspace_encoding_step_2" => Limit(0, Nz-1, ceil(Int, Nz / 2)-1),#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, Ns-1, ceil(Int, Ns / 2)-1), #min, max, center, e.g. imaging slice number + "enc_lim_slice" => Limit(0, Ns-1, ceil(Int, Ns / 2)-1),#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 diff --git a/src/KomaUI.jl b/src/KomaUI.jl index fb52033cd..d5c5d38b9 100644 --- a/src/KomaUI.jl +++ b/src/KomaUI.jl @@ -217,19 +217,24 @@ function KomaUI(; darkmode=true, frame=true, phantom_mode="2D", sim=Dict{String, # Get the value of the raw signal and prepare for reconstruction raw_aux = raw_ui[] + traj_dims = convert(Int32, raw_aux.profiles[1].head.trajectory_dimensions) raw_aux.profiles = raw_aux.profiles[getproperty.(getproperty.(raw_aux.profiles, :head), :flags) .!= 268435456] #Extra profile in JEMRIS simulations acq_data = AcquisitionData(raw_aux) acq_data.traj[1].circular = false #Removing circular window acq_data.traj[1].nodes = acq_data.traj[1].nodes[1:2,:] ./ maximum(2*abs.(acq_data.traj[1].nodes[:])) #Normalize k-space to -.5 to .5 for NUFFT - Nx, Ny, Nz = raw_aux.params["reconSize"][1:3] - rec_params[:reconSize] = (Nx, Ny, Nz) - rec_params[:densityWeighting] = true + Nx, Ny, Nz = raw_aux.params["reconSize"] # Perform reconstruction - @info "Running reconstruction ..." + rec_params[:densityWeighting] = true + if traj_dims == 2 + rec_params[:reconSize] = (Nx, Ny) + @info "Running 2D reconstruction ..." + else + rec_params[:reconSize] = (Nx, Ny, Nz) + @info "Running 2D reconstruction ..." + end rec_aux = @timed reconstruction(acq_data, rec_params) image = reshape(rec_aux.value.data, Nx, Ny, :) - # After Recon go to Image recon_time = rec_aux.time @js_ w (document.getElementById("recon!").innerHTML = "Reconstruct!") From 570bde79d30451b966e9fcecdbd752cf54c220dd Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Sat, 8 Jun 2024 13:51:52 -0500 Subject: [PATCH 13/30] Debugging reconstruction calls from KomaUI(). modified: KomaMRICore/src/rawdata/ISMRMRD.jl modified: src/KomaUI.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 6 +++++- src/KomaUI.jl | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index bf7b2983b..56a2712bc 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -111,7 +111,11 @@ function signal_to_raw_data( else ndims = Nd_seq end - @info "Creating $ndims dimensional ISMRMRD file ..." + if ndims == 2 + @info "Creating 2D ISMRMRD file ..." + else + @info "Creating 3D 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[:])) diff --git a/src/KomaUI.jl b/src/KomaUI.jl index d5c5d38b9..ca87d8693 100644 --- a/src/KomaUI.jl +++ b/src/KomaUI.jl @@ -231,7 +231,7 @@ function KomaUI(; darkmode=true, frame=true, phantom_mode="2D", sim=Dict{String, @info "Running 2D reconstruction ..." else rec_params[:reconSize] = (Nx, Ny, Nz) - @info "Running 2D reconstruction ..." + @info "Running 3D reconstruction ..." end rec_aux = @timed reconstruction(acq_data, rec_params) image = reshape(rec_aux.value.data, Nx, Ny, :) From 3bd4209daf7552daa1eb0d0ffd41d40a24e59650 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Sat, 8 Jun 2024 16:11:55 -0500 Subject: [PATCH 14/30] Debug raw data conversion for ge3d_cac.seq, something funny...!!! modified: KomaMRICore/src/rawdata/ISMRMRD.jl modified: src/KomaUI.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 24 ++++++++++++++---------- src/KomaUI.jl | 13 ++++++++----- 2 files changed, 22 insertions(+), 15 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 56a2712bc..281003efb 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -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) Transforms the raw signal into a RawAcquisitionData struct (nearly equivalent to the ISMRMRD format) used for reconstruction with MRIReco. @@ -56,6 +56,7 @@ 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 # Returns - `raw`: (`::RawAcquisitionData`) RawAcquisitionData struct @@ -70,15 +71,18 @@ 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=-1, + phantom_name="Phantom", sys=Scanner(), sim_params=Dict{String,Any}(), ndims=0, ) + + @debug "Something going funny here, making 960 extra profiles for ge3d_cac.seq! *** CAC 240608" + version = string(VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "..", "Project.toml"))["version"])) #Number of samples and FOV _, ktraj = get_kspace(seq) #kspace information @@ -106,7 +110,7 @@ function signal_to_raw_data( FOVz = Ny * Δx[3] end Nd_seq = (Nx > 1) + (Ny > 1) + (Nz > 1) - if ndims > -1 #just in case ndims is provided in a call + if ndims < 0 #just in case ndims is provided in a call if ndims != Nd_seq; @warn("Seqfile is $Nd_seq dimensional but .mrd file will be $ndims."); end else ndims = Nd_seq @@ -148,11 +152,11 @@ function signal_to_raw_data( "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)-1),#min, max, center, e.g. phase encoding line number - "enc_lim_kspace_encoding_step_1" => Limit(0, Ny-1, ceil(Int, Ny / 2)-1),#min, max, center, e.g. partition encoding number - "enc_lim_kspace_encoding_step_2" => Limit(0, Nz-1, ceil(Int, Nz / 2)-1),#min, max, center, e.g. partition encoding number + "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, Nz-1, ceil(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, Ns-1, ceil(Int, Ns / 2)-1),#min, max, center, e.g. imaging slice number + "enc_lim_slice" => Limit(0, Ns-1, ceil(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 @@ -240,11 +244,11 @@ function signal_to_raw_data( current += Nsamples if scan_counter % NadcsPerSlice == 0 ns += 1 #another slice - #scan_counter = 0 #reset counter + scan_counter = 0 #reset counter end if scan_counter % NadcsPerPE1 == 0 #Using Ns for slice number nz += 1 #another kspace_encode_step_2 - #scan_counter = 0 #reset counter + scan_counter = 0 #reset counter end # more here for other counters, i.e. dynamic end diff --git a/src/KomaUI.jl b/src/KomaUI.jl index ca87d8693..bc6d09cc1 100644 --- a/src/KomaUI.jl +++ b/src/KomaUI.jl @@ -219,20 +219,23 @@ function KomaUI(; darkmode=true, frame=true, phantom_mode="2D", sim=Dict{String, raw_aux = raw_ui[] traj_dims = convert(Int32, raw_aux.profiles[1].head.trajectory_dimensions) raw_aux.profiles = raw_aux.profiles[getproperty.(getproperty.(raw_aux.profiles, :head), :flags) .!= 268435456] #Extra profile in JEMRIS simulations + Nx, Ny, Nz = raw_aux.params["reconSize"] acq_data = AcquisitionData(raw_aux) acq_data.traj[1].circular = false #Removing circular window - acq_data.traj[1].nodes = acq_data.traj[1].nodes[1:2,:] ./ maximum(2*abs.(acq_data.traj[1].nodes[:])) #Normalize k-space to -.5 to .5 for NUFFT - Nx, Ny, Nz = raw_aux.params["reconSize"] - - # Perform reconstruction rec_params[:densityWeighting] = true if traj_dims == 2 + acq_data.traj[1].nodes = acq_data.traj[1].nodes[1:2,:] ./ maximum(2*abs.(acq_data.traj[1].nodes[:])) #Normalize k-space to -.5 to .5 for NUFFT rec_params[:reconSize] = (Nx, Ny) @info "Running 2D reconstruction ..." - else + else + acq_data.traj[1].nodes = acq_data.traj[1].nodes ./ maximum(2*abs.(acq_data.traj[1].nodes[:])) #Normalize k-space to -.5 to .5 for NUFFT rec_params[:reconSize] = (Nx, Ny, Nz) @info "Running 3D reconstruction ..." end + nodes_size = size( acq_data.traj[1].nodes) + @debug rec_params nodes_size + + # Perform reconstruction rec_aux = @timed reconstruction(acq_data, rec_params) image = reshape(rec_aux.value.data, Nx, Ny, :) # After Recon go to Image From d977f2529d4655b6db14978a985659bee9b5745a Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Sun, 9 Jun 2024 22:22:53 -0500 Subject: [PATCH 15/30] A fair amount of work to estimate proper 3D recon parameters from seq, and work and debugging of the scan counters for 3D. modified: KomaMRICore/src/rawdata/ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 172 +++++++++++++++++++---------- 1 file changed, 112 insertions(+), 60 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 281003efb..30f49b464 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -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, ndims) + raw = signal_to_raw_data(signal, seq; phantom_name, sys, sim_params, ndims=2) Transforms the raw signal into a RawAcquisitionData struct (nearly equivalent to the ISMRMRD format) used for reconstruction with MRIReco. @@ -78,48 +78,25 @@ julia> plot_signal(raw) """ function signal_to_raw_data( signal, seq; - phantom_name="Phantom", sys=Scanner(), sim_params=Dict{String,Any}(), ndims=0, + phantom_name="Phantom", sys=Scanner(), sim_params=Dict{String,Any}(), ndims=2, use_ndseq=true ) - - @debug "Something going funny here, making 960 extra profiles for ge3d_cac.seq! *** CAC 240608" - 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 - idxs_zero = findall(iszero, Wk) #check for zeros - Wk[idxs_zero] .= 1e6 - Δx = 1 ./ Wk[1:3] #[m] x-y-z - Nx = get(seq.DEF, "Nx", 1) - Ny = get(seq.DEF, "Ny", 1) - Nz = get(seq.DEF, "Nz", 1) - Ns = get(seq.DEF, "Ns", 1) # number of slices or slabs - if haskey(seq.DEF, "FOV") - FOVx, FOVy, FOVz = 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 - if FOVz > 1 FOVz *= 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]) - Nz = round(Int64, FOVy / Δx[3]) + + #Estimate sequence reconstruction dimension, using seq.DEF and potentially other information + Nd_seq, Nx, Ny, Nz, Ns, FOVx, FOVy, FOVz, Δx, ktraj = estimate_seq_recon_dimension(signal, 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 - FOVx = Nx * Δx[1] - FOVy = Ny * Δx[2] - FOVz = Ny * Δx[3] - end - Nd_seq = (Nx > 1) + (Ny > 1) + (Nz > 1) - if ndims < 0 #just in case ndims is provided in a call if ndims != Nd_seq; @warn("Seqfile is $Nd_seq dimensional but .mrd file will be $ndims."); end - else - ndims = Nd_seq end if ndims == 2 @info "Creating 2D ISMRMRD file ..." - else + elseif ndims == 3 @info "Creating 3D ISMRMRD file ..." - end + else + @warn("Seqfile is $Nd_seq dimensional which is not yet supported.") + end + #It needs to be transposed for the raw data ktraj = maximum(2*abs.(ktraj[:])) == 0 ? transpose(ktraj) : transpose(ktraj)./ maximum(2*abs.(ktraj[:])) @@ -152,11 +129,11 @@ function signal_to_raw_data( "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, Nz-1, ceil(Int, Nz / 2)),#min, max, center, e.g. partition encoding number + "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, Nz-1, ceil(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, Ns-1, ceil(Int, Ns / 2)),#min, max, center, e.g. imaging slice number + "enc_lim_slice" => Limit(0, Ns-1, ceil(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 @@ -175,13 +152,15 @@ 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)) - NadcsPerSlice = floor(Int, Nadcs / Ns) +1 - NadcsPerPE1 = floor(Int, Nadcs / Nz) +1 + Nro = sum(is_ADC_on.(seq)) + NroPerPE1 = floor(Int, Nro / (Ny*Ns*Nz)) + NroPerSlice = floor(Int, Nro / Ns) + NroPerPE2 = floor(Int, Nro / Nz) scan_counter = 0 - ns = 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] @@ -191,9 +170,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 - flag += ISMRMRD_ACQ_LAST_IN_ENCODE_STEP1 - flag += ISMRMRD_ACQ_LAST_IN_SLICE + 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( @@ -219,15 +200,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(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(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 + 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 @@ -240,15 +221,17 @@ 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 % NadcsPerSlice == 0 + if scan_counter % NroPerSlice == 0 ns += 1 #another slice - scan_counter = 0 #reset counter end - if scan_counter % NadcsPerPE1 == 0 #Using Ns for slice number + if scan_counter % NroPerPE1 == 0 + ny += 1 #another kspace_encode_step_2 + end + if scan_counter % NroPerPE2 == 0 + ny = 0 nz += 1 #another kspace_encode_step_2 - scan_counter = 0 #reset counter end # more here for other counters, i.e. dynamic end @@ -257,6 +240,75 @@ function signal_to_raw_data( return RawAcquisitionData(params, profiles) end +""" + Nd_seq, Nx, Ny, Nz, Ns, FOVx, FOVy, FOVz, Δx = estimate_seq_recon_dimension(signal, seq; sim_params) + +Utility function for the best estimate of the reconstruction dimension. + +# Arguments +- `signal`: (`::Matrix{Complex}`) raw signal matrix +- `seq`: (`::Sequence`) Sequence struct + +# Keywords +- `sim_params`: (`::Dict{String, Any}`, `=Dict{String,Any}()`) simulation parameter dictionary + +# 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> signal = simulate(obj, seq, sys; sim_params) + +julia> Nd_seq = estimate_seq_recon_dimension(seq, signal; sim_params) + +julia> OUTPUT **** +``` +""" +function estimate_seq_recon_dimension( + signal, seq; sim_params=Dict{String,Any}(), +) + + #Number of samples and FOV + _, 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[idxs_zero] .= 1e6 + Δx = 1 ./ Wk[1:3] #[m] x-y-z + 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 haskey(seq.DEF, "FOV") + FOVx, FOVy, FOVz = 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 + if FOVz > 1 FOVz *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm + 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, FOVy / Δx[3]) end + else + @warn "Estimating FOV parameters from seq.DEF Nx, Ny, Nz and traj." + FOVx = Nx * Δx[1] + FOVy = Ny * Δx[2] + FOVz = Ny * Δx[3] + end + Nd_seq = (Nx > 1) + (Ny > 1) + (Nz > 1) + 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) From a4c387c92612e4962015255ad273b352971e0ebc Mon Sep 17 00:00:00 2001 From: Curt Corum <1796102+curtcorum@users.noreply.github.com> Date: Wed, 12 Jun 2024 16:05:45 -0500 Subject: [PATCH 16/30] Update ISMRMRD.jl Should fix test failure at: https://github.com/JuliaHealth/KomaMRI.jl/blob/931fec8f318b247dfbc90fa33f0b598216f6a3c8/KomaMRICore/test/runtests.jl#L154-L158 Mentioned: https://github.com/JuliaHealth/KomaMRI.jl/pull/324#issuecomment-2160854813 In draft PR: 3D recon #324 --- KomaMRICore/src/rawdata/ISMRMRD.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 30f49b464..9f763527a 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -78,7 +78,7 @@ julia> plot_signal(raw) """ function signal_to_raw_data( signal, seq; - phantom_name="Phantom", sys=Scanner(), sim_params=Dict{String,Any}(), ndims=2, use_ndseq=true + phantom_name="Phantom", sys=Scanner(), sim_params=Dict{String,Any}(), ndims=2, use_ndseq=false ) version = string(VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "..", "Project.toml"))["version"])) From bd51a0166fe5d5d2c2667e0b62999fb435545214 Mon Sep 17 00:00:00 2001 From: Curt Corum <1796102+curtcorum@users.noreply.github.com> Date: Wed, 12 Jun 2024 16:44:09 -0500 Subject: [PATCH 17/30] Update ISMRMRD.jl Added `use_ndseq` to `signal_to_raw_data()` docstring --- KomaMRICore/src/rawdata/ISMRMRD.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 9f763527a..1a03a42c7 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -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, ndims=2) + 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. @@ -57,6 +57,7 @@ format) used for reconstruction with MRIReco. - `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 From 9aac620e60c499d96eacb72818ca80aec3098e9a Mon Sep 17 00:00:00 2001 From: Curt Corum <1796102+curtcorum@users.noreply.github.com> Date: Wed, 12 Jun 2024 16:57:59 -0500 Subject: [PATCH 18/30] Update ISMRMRD.jl More updates for FOV estimation --- KomaMRICore/src/rawdata/ISMRMRD.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 1a03a42c7..60862dcfb 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -289,17 +289,17 @@ function estimate_seq_recon_dimension( Ns = ceil(Int64, get(seq.DEF, "Ns", 1)) # number of slices or slabs if haskey(seq.DEF, "FOV") FOVx, FOVy, FOVz = 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 - if FOVz > 1 FOVz *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm - 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, FOVy / Δx[3]) end + 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 + if FOVz > 1 FOVz *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm + if Nx < 2 Nx = ceil(Int64, FOVx / Δx[1]) end + if Ny < 2 Ny = ceil(Int64, FOVy / Δx[2]) end + if Nz < 2 Nz = ceil(Int64, FOVy / Δx[3]) end else @warn "Estimating FOV parameters from seq.DEF Nx, Ny, Nz and traj." - FOVx = Nx * Δx[1] - FOVy = Ny * Δx[2] - FOVz = Ny * Δx[3] + if Nx > 1 FOVx = Nx * Δx[1] else FOVx = 0 end + if Ny > 1 FOVy = Ny * Δx[2] else FOVy = 0 end + if Ny > 1 FOVz = Nz * Δx[3] else FOVz = 0 end end Nd_seq = (Nx > 1) + (Ny > 1) + (Nz > 1) s_ktraj = size( ktraj) From 05d862a69031bd03c9030723ea93bd1cecbced0c Mon Sep 17 00:00:00 2001 From: Curt Corum <1796102+curtcorum@users.noreply.github.com> Date: Wed, 12 Jun 2024 17:09:03 -0500 Subject: [PATCH 19/30] Update ISMRMRD.jl debugging --- KomaMRICore/src/rawdata/ISMRMRD.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 60862dcfb..2f54eba05 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -292,9 +292,9 @@ function estimate_seq_recon_dimension( 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 if FOVz > 1 FOVz *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm - if Nx < 2 Nx = ceil(Int64, FOVx / Δx[1]) end - if Ny < 2 Ny = ceil(Int64, FOVy / Δx[2]) end - if Nz < 2 Nz = ceil(Int64, FOVy / Δx[3]) end + 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, FOVy / Δx[3]) end else @warn "Estimating FOV parameters from seq.DEF Nx, Ny, Nz and traj." if Nx > 1 FOVx = Nx * Δx[1] else FOVx = 0 end From 8719f0a93f992f2e2503459161796f57294e919b Mon Sep 17 00:00:00 2001 From: Curt Corum <1796102+curtcorum@users.noreply.github.com> Date: Wed, 12 Jun 2024 17:12:38 -0500 Subject: [PATCH 20/30] Update ISMRMRD.jl One more debug? --- KomaMRICore/src/rawdata/ISMRMRD.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 2f54eba05..ddb7a280e 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -297,9 +297,9 @@ function estimate_seq_recon_dimension( if Nz == 1 Nz = ceil(Int64, FOVy / Δx[3]) end else @warn "Estimating FOV parameters from seq.DEF Nx, Ny, Nz and traj." - if Nx > 1 FOVx = Nx * Δx[1] else FOVx = 0 end - if Ny > 1 FOVy = Ny * Δx[2] else FOVy = 0 end - if Ny > 1 FOVz = Nz * Δx[3] else FOVz = 0 end + if Nx >= 1 FOVx = Nx * Δx[1] else FOVx = 0 end + if Ny >= 1 FOVy = Ny * Δx[2] else FOVy = 0 end + if Ny >= 1 FOVz = Nz * Δx[3] else FOVz = 0 end end Nd_seq = (Nx > 1) + (Ny > 1) + (Nz > 1) s_ktraj = size( ktraj) From acc981dc49f9410c29086f0d579b1007c7b6a395 Mon Sep 17 00:00:00 2001 From: Curt Corum <1796102+curtcorum@users.noreply.github.com> Date: Wed, 12 Jun 2024 17:27:45 -0500 Subject: [PATCH 21/30] Update ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index ddb7a280e..4e3fcb5a8 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -95,7 +95,7 @@ function signal_to_raw_data( elseif ndims == 3 @info "Creating 3D ISMRMRD file ..." else - @warn("Seqfile is $Nd_seq dimensional which is not yet supported.") + @info("Creating a $ndims dimensional ISMRMRD file...") end #It needs to be transposed for the raw data @@ -281,7 +281,7 @@ function estimate_seq_recon_dimension( maxk = maximum(ktraj, dims=1) Wk = maxk .- mink idxs_zero = findall(iszero, Wk) #check for zeros - Wk[idxs_zero] .= 1e6 + Wk[idxs_zero] .= 1.0e6 Δx = 1 ./ Wk[1:3] #[m] x-y-z Nx = ceil(Int64, get(seq.DEF, "Nx", 1)) Ny = ceil(Int64, get(seq.DEF, "Ny", 1)) @@ -292,14 +292,14 @@ function estimate_seq_recon_dimension( 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 if FOVz > 1 FOVz *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm - 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, FOVy / Δx[3]) end + #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, FOVy / Δx[3]) end else @warn "Estimating FOV parameters from seq.DEF Nx, Ny, Nz and traj." - if Nx >= 1 FOVx = Nx * Δx[1] else FOVx = 0 end - if Ny >= 1 FOVy = Ny * Δx[2] else FOVy = 0 end - if Ny >= 1 FOVz = Nz * Δx[3] else FOVz = 0 end + if ~idx_zero[1] FOVx = Nx * Δx[1] else FOVx = 0 end + if ~idx_zero[2] FOVy = Ny * Δx[2] else FOVy = 0 end + if ~idx_zero[3] FOVz = Nz * Δx[3] else FOVz = 0 end end Nd_seq = (Nx > 1) + (Ny > 1) + (Nz > 1) s_ktraj = size( ktraj) From eb9b789656f58f4992730b663d20dc2ea3075dcf Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Wed, 12 Jun 2024 18:05:53 -0500 Subject: [PATCH 22/30] Debugging signal_to_raw_data modified: KomaMRICore/src/rawdata/ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 4e3fcb5a8..fc246d6e0 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -281,25 +281,26 @@ function estimate_seq_recon_dimension( 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.0e6 Δx = 1 ./ Wk[1:3] #[m] x-y-z 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 haskey(seq.DEF, "FOV") + if haskey(seq.DEF, "FOV") #needs info or warning for other substitutions??? FOVx, FOVy, FOVz = 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 if FOVz > 1 FOVz *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm - #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, FOVy / Δx[3]) end + 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 else @warn "Estimating FOV parameters from seq.DEF Nx, Ny, Nz and traj." - if ~idx_zero[1] FOVx = Nx * Δx[1] else FOVx = 0 end - if ~idx_zero[2] FOVy = Ny * Δx[2] else FOVy = 0 end - if ~idx_zero[3] FOVz = Nz * Δx[3] else FOVz = 0 end + if !Wk_el_iszero[1] FOVx = Nx * Δx[1] else FOVx = 0 end + if !Wk_el_iszero[2] FOVy = Ny * Δx[2] else FOVy = 0 end + if !Wk_el_iszero[3] FOVz = Nz * Δx[3] else FOVz = 0 end end Nd_seq = (Nx > 1) + (Ny > 1) + (Nz > 1) s_ktraj = size( ktraj) From 9f25e2499d2aa041d5cc144f11db1e04d7aeacff Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Wed, 12 Jun 2024 18:35:57 -0500 Subject: [PATCH 23/30] Check on Ny, Ns, Nz == 0 for scan counters. modified: KomaMRICore/src/rawdata/ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index fc246d6e0..0c54b656a 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -154,9 +154,10 @@ function signal_to_raw_data( profiles = Profile[] t_acq = get_adc_sampling_times(seq) Nro = sum(is_ADC_on.(seq)) - NroPerPE1 = floor(Int, Nro / (Ny*Ns*Nz)) - NroPerSlice = floor(Int, Nro / Ns) - NroPerPE2 = floor(Int, Nro / Nz) + Nyt = max(Ny, 1); Nst = max(Ns, 1); Nzt = max(Nz, 1) #zero is really one... + NroPerPE1 = floor(Int, Nro / (Nyt*Nst*Nzt)) + NroPerSlice = floor(Int, Nro / Nst) + NroPerPE2 = floor(Int, Nro / Nzt) scan_counter = 0 ny = 0 #PE1 counter ns = 0 #Slice counter @@ -298,9 +299,9 @@ function estimate_seq_recon_dimension( if Nz == 1 Nz = ceil(Int64, FOVz / Δx[3]) end else @warn "Estimating FOV parameters from seq.DEF Nx, Ny, Nz and traj." - if !Wk_el_iszero[1] FOVx = Nx * Δx[1] else FOVx = 0 end - if !Wk_el_iszero[2] FOVy = Ny * Δx[2] else FOVy = 0 end - if !Wk_el_iszero[3] FOVz = Nz * Δx[3] else FOVz = 0 end + if !Wk_el_iszero[1] FOVx = Nx * Δx[1] else FOVx = 1 end + if !Wk_el_iszero[2] FOVy = Ny * Δx[2] else FOVy = 1 end + if !Wk_el_iszero[3] FOVz = Nz * Δx[3] else FOVz = 1 end end Nd_seq = (Nx > 1) + (Ny > 1) + (Nz > 1) s_ktraj = size( ktraj) From 26c6c1b9140b2fb7c47cf1b8aec4f186cb61a47f Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Thu, 13 Jun 2024 16:45:46 -0500 Subject: [PATCH 24/30] Defaulting to `use_ndseq=true` for `signal_to_raw_data`. modified: KomaMRICore/src/rawdata/ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 2 +- KomaMRICore/src/simulation/GPUFunctions.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 0c54b656a..23a7d20c2 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -79,7 +79,7 @@ julia> plot_signal(raw) """ function signal_to_raw_data( signal, seq; - phantom_name="Phantom", sys=Scanner(), sim_params=Dict{String,Any}(), ndims=2, use_ndseq=false + phantom_name="Phantom", sys=Scanner(), sim_params=Dict{String,Any}(), ndims=2, use_ndseq=true ) version = string(VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "..", "Project.toml"))["version"])) diff --git a/KomaMRICore/src/simulation/GPUFunctions.jl b/KomaMRICore/src/simulation/GPUFunctions.jl index 6d6200678..607509e83 100644 --- a/KomaMRICore/src/simulation/GPUFunctions.jl +++ b/KomaMRICore/src/simulation/GPUFunctions.jl @@ -92,4 +92,4 @@ function print_devices(use_gpu = true) _print_devices(backend) end -export print_devices \ No newline at end of file +export print_devices From 585304f6d95e543837b11ade7efe8cccf31ed333 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Mon, 8 Jul 2024 22:05:16 -0500 Subject: [PATCH 25/30] Working on central recon parameter estimation in KomaMRICore/src/rawdata/ISMRMRD.jl modified: KomaMRICore/src/KomaMRICore.jl modified: KomaMRICore/src/rawdata/ISMRMRD.jl modified: KomaMRIFiles/src/Sequence/Pulseq.jl --- KomaMRICore/src/KomaMRICore.jl | 3 ++- KomaMRICore/src/rawdata/ISMRMRD.jl | 33 +++++++++++++++++++++++------ KomaMRIFiles/src/Sequence/Pulseq.jl | 8 ++++--- 3 files changed, 33 insertions(+), 11 deletions(-) diff --git a/KomaMRICore/src/KomaMRICore.jl b/KomaMRICore/src/KomaMRICore.jl index 8f3e33333..ccbf68468 100644 --- a/KomaMRICore/src/KomaMRICore.jl +++ b/KomaMRICore/src/KomaMRICore.jl @@ -22,7 +22,8 @@ include("simulation/Functors.jl") include("simulation/SimulatorCore.jl") # ISMRMRD -export signal_to_raw_data +#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, default_sim_params diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 23a7d20c2..1ffd3c8b7 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -84,7 +84,7 @@ function signal_to_raw_data( version = string(VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "..", "Project.toml"))["version"])) #Estimate sequence reconstruction dimension, using seq.DEF and potentially other information - Nd_seq, Nx, Ny, Nz, Ns, FOVx, FOVy, FOVz, Δx, ktraj = estimate_seq_recon_dimension(signal, seq; sim_params) + 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 @@ -243,12 +243,11 @@ function signal_to_raw_data( end """ - Nd_seq, Nx, Ny, Nz, Ns, FOVx, FOVy, FOVz, Δx = estimate_seq_recon_dimension(signal, seq; sim_params) + Nd_seq, Nx, Ny, Nz, Ns, FOVx, FOVy, FOVz, Δx = estimate_seq_recon_dimension(seq; sim_params) Utility function for the best estimate of the reconstruction dimension. # Arguments -- `signal`: (`::Matrix{Complex}`) raw signal matrix - `seq`: (`::Sequence`) Sequence struct # Keywords @@ -272,11 +271,11 @@ julia> Nd_seq = estimate_seq_recon_dimension(seq, signal; sim_params) julia> OUTPUT **** ``` """ -function estimate_seq_recon_dimension( - signal, seq; sim_params=Dict{String,Any}(), +function estimate_seq_recon_dimension(seq; sim_params=Dict{String,Any}(), ) - - #Number of samples and FOV + #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) @@ -285,6 +284,26 @@ function estimate_seq_recon_dimension( Wk_el_iszero = iszero.( Wk) # bool array for later Wk[idxs_zero] .= 1.0e6 Δx = 1 ./ Wk[1:3] #[m] x-y-z + + + #estimate sequence acquisition structure + Ns_seq = length(unique(seq.RF.Δf)) #slices, slabs, or Cartesean off-center, 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" + + # Guessing Cartesean recon dimensions + # ideally all estimates of recon dimensions in one place, as late as possible, 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 + + # some guesses + seq_3d=false; seq_2d=false; seq_nd=false; seq_cartesean=false; seq_radial=false; seq_spiral=false + if Ns_seq == Nv_seq; seq_rad = true; 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)) diff --git a/KomaMRIFiles/src/Sequence/Pulseq.jl b/KomaMRIFiles/src/Sequence/Pulseq.jl index ef67baefc..25f8e8c6b 100644 --- a/KomaMRIFiles/src/Sequence/Pulseq.jl +++ b/KomaMRIFiles/src/Sequence/Pulseq.jl @@ -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 From af870bdcd83f102bef3175c9c4dd4c405ea6c0ff Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Tue, 16 Jul 2024 00:03:22 -0500 Subject: [PATCH 26/30] Working on estimate_seq_recon_dimension, working for 2d multislice, not for spiral modified: KomaMRICore/src/rawdata/ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 138 ++++++++++++++++++++--------- 1 file changed, 94 insertions(+), 44 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 1ffd3c8b7..c60e17fa9 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -130,11 +130,11 @@ function signal_to_raw_data( "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, Nz-1, ceil(Int, Nz / 2)), #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, Ns-1, ceil(Int, Ns / 2)), #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 @@ -154,10 +154,11 @@ function signal_to_raw_data( profiles = Profile[] t_acq = get_adc_sampling_times(seq) Nro = sum(is_ADC_on.(seq)) - Nyt = max(Ny, 1); Nst = max(Ns, 1); Nzt = max(Nz, 1) #zero is really one... - NroPerPE1 = floor(Int, Nro / (Nyt*Nst*Nzt)) - NroPerSlice = floor(Int, Nro / Nst) - NroPerPE2 = floor(Int, Nro / Nzt) + #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 ny = 0 #PE1 counter ns = 0 #Slice counter @@ -174,9 +175,9 @@ function signal_to_raw_data( flag += ISMRMRD_ACQ_FIRST_IN_SLICE 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 + 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( @@ -225,20 +226,22 @@ function signal_to_raw_data( #Update counters scan_counter += 1 #another ro current += Nsamples - if scan_counter % NroPerSlice == 0 - ns += 1 #another slice - end if scan_counter % NroPerPE1 == 0 - ny += 1 #another kspace_encode_step_2 + ny += 1 #another kspace_encode_step_1 end + ny = ny % Ny if scan_counter % NroPerPE2 == 0 - ny = 0 nz += 1 #another kspace_encode_step_2 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 @@ -251,7 +254,7 @@ Utility function for the best estimate of the reconstruction dimension. - `seq`: (`::Sequence`) Sequence struct # Keywords -- `sim_params`: (`::Dict{String, Any}`, `=Dict{String,Any}()`) simulation parameter dictionary +- `sim_params`: (`::Dict{String, Any}`, `=Dict{String,Any}()`) simulation parameter dictionary (IGNORED) # Returns - `Nd_seq`: (`Int`) Estimated reconstruction dimension of seq. @@ -282,47 +285,94 @@ function estimate_seq_recon_dimension(seq; sim_params=Dict{String,Any}(), Wk = maxk .- mink idxs_zero = findall(iszero, Wk) #check for zeros Wk_el_iszero = iszero.( Wk) # bool array for later - Wk[idxs_zero] .= 1.0e6 + Wk[idxs_zero] .= 1.0e-6 Δx = 1 ./ Wk[1:3] #[m] x-y-z - #estimate sequence acquisition structure Ns_seq = length(unique(seq.RF.Δf)) #slices, slabs, or Cartesean off-center, 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" - # Guessing Cartesean recon dimensions - # ideally all estimates of recon dimensions in one place, as late as possible, 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 + #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_3d=false; seq_2d=false; seq_nd=false; seq_cartesean=false; seq_radial=false; seq_spiral=false - if Ns_seq == Nv_seq; seq_rad = true; end + 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 Ns_seq == Nv_seq + 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 = get(seq.DEF, "Ns", Ns_seq) #slices or slabs + Nx = get(seq.DEF, "Nx", Np_seq) + if seq_cartesean + if seq_2d + Ny = get(seq.DEF, "Ny", ceil(Int64, Nv_seq/Ns)) #pe1 + Nz = 1 + elseif seq_3d + Ny = get(seq.DEF, "Ny", Ny_k) #pe1 + Nz = get(seq.DEF, "Nz", ceil(Int64, Nv_seq/Ny_k)) #pe2 + end + else + Ny = ceil(Int64, Nv_seq/Ns) + 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 + #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 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 - if FOVz > 1 FOVz *= 1e-3 end #mm to m, older versions of Pulseq saved FOV in mm - 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 - else - @warn "Estimating FOV parameters from seq.DEF Nx, Ny, Nz and traj." - if !Wk_el_iszero[1] FOVx = Nx * Δx[1] else FOVx = 1 end - if !Wk_el_iszero[2] FOVy = Ny * Δx[2] else FOVy = 1 end - if !Wk_el_iszero[3] FOVz = Nz * Δx[3] else FOVz = 1 end + if FOVx > 1.0 + FOVx *= 1e-3 + FOVy *= 1e-3 + FOVz *= 1e-3 + @warn "Scaling FOV to m from older pulseq file mm." + 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 = (Nx > 1) + (Ny > 1) + (Nz > 1) + 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" From 12e0d63e64e904ee026601caab49c45c6680edda Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Tue, 16 Jul 2024 20:06:46 -0500 Subject: [PATCH 27/30] Fixing merge. modified: KomaMRICore/src/rawdata/ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index fcf5e551c..4bd435a04 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -81,7 +81,7 @@ function signal_to_raw_data( signal, seq; phantom_name="Phantom", sys=Scanner(), sim_params=Dict{String,Any}(), ndims=2, use_ndseq=true ) - version = string(VersionNumber(Pkg.TOML.parsefile(joinpath(@__DIR__, "..", "..", "Project.toml"))["version"])) + #Number of samples and FOV Nd_seq, Nx, Ny, Nz, Ns, FOVx, FOVy, FOVz, Δx, ktraj = estimate_seq_recon_dimension(seq; sim_params) if use_ndseq @@ -111,7 +111,7 @@ function signal_to_raw_data( params = Dict( #AcquisitionSystemInformation "systemVendor" => "KomaMRI.jl", #String - "systemModel" => "v"*version, #String + "systemModel" => string(pkgversion(@__MODULE__)), #String "systemFieldStrength_T" => sys.B0, #Float "institutionName" => "Pontificia Universidad Catolica de Chile", #String #subjectInformation @@ -203,7 +203,7 @@ function signal_to_raw_data( Float32.((0, 0, 0)), #patient_table_position float32x3: Patient table off-center EncodingCounters( #idx uint16x17: Encoding loop counters UInt16(ny), #kspace_encode_step_1 uint16: e.g. phase encoding line number - UInt16(nz), #slice uint16: e.g. imaging slice 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 From b101fe1700787ccc1dabbbfb80548e723144ccc3 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Tue, 16 Jul 2024 20:34:41 -0500 Subject: [PATCH 28/30] Whitespace formatting. modified: KomaMRICore/src/rawdata/ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 4bd435a04..633dd8e2f 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -129,11 +129,11 @@ function signal_to_raw_data( "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, 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_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, Ns-1, floor(Int, Ns / 2)), #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 From c98d92b1de2bb170d04aa3c1ba9626de2f5d6665 Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Tue, 16 Jul 2024 21:30:43 -0500 Subject: [PATCH 29/30] Forced type of Nx, Ny, Nz, Ns read from seq.DEF to be Int64. Was not matching ::Tuple{Vararg{Int64, D}} in MRIBase ~/.julia/packages/MRIBase/GskMD/src/Datatypes/AcqData.jl:21 modified: KomaMRICore/src/rawdata/ISMRMRD.jl --- KomaMRICore/src/rawdata/ISMRMRD.jl | 18 +++++++----------- 1 file changed, 7 insertions(+), 11 deletions(-) diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 633dd8e2f..45b4e9ea7 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -245,7 +245,7 @@ function signal_to_raw_data( end """ - Nd_seq, Nx, Ny, Nz, Ns, FOVx, FOVy, FOVz, Δx = estimate_seq_recon_dimension(seq; sim_params) + 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. @@ -266,11 +266,7 @@ julia> sys, obj, seq = Scanner(), brain_phantom2D(), read_seq(seq_file) julia> sim_params = KomaMRICore.default_sim_params(); sim_params["return_type"] = "mat" -julia> signal = simulate(obj, seq, sys; sim_params) - -julia> Nd_seq = estimate_seq_recon_dimension(seq, signal; sim_params) - -julia> OUTPUT **** +julia> Nd_seq = estimate_seq_recon_dimension(seq; sim_params) ``` """ function estimate_seq_recon_dimension(seq; sim_params=Dict{String,Any}(), @@ -323,15 +319,15 @@ function estimate_seq_recon_dimension(seq; sim_params=Dict{String,Any}(), # Guessing Cartesean recon dimensions # ideally all estimates of recon dimensions in one place, as late as possible, non-Cartesean ?? *** CAC 240708 - Ns = get(seq.DEF, "Ns", Ns_seq) #slices or slabs - Nx = get(seq.DEF, "Nx", Np_seq) + 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 = get(seq.DEF, "Ny", ceil(Int64, Nv_seq/Ns)) #pe1 + Ny = Int64(get(seq.DEF, "Ny", ceil(Int64, Nv_seq/Ns))) #pe1 Nz = 1 elseif seq_3d - Ny = get(seq.DEF, "Ny", Ny_k) #pe1 - Nz = get(seq.DEF, "Nz", ceil(Int64, Nv_seq/Ny_k)) #pe2 + Ny = Int64(get(seq.DEF, "Ny", Ny_k)) #pe1 + Nz = Int64(get(seq.DEF, "Nz", ceil(Int64, Nv_seq/Ny_k))) #pe2 end else Ny = ceil(Int64, Nv_seq/Ns) From 0f317cb3eae59f59747a00ee8cc85de3420323eb Mon Sep 17 00:00:00 2001 From: curtcorum <1796102+curtcorum@users.noreply.github.com> Date: Wed, 17 Jul 2024 21:26:37 -0500 Subject: [PATCH 30/30] LiearAlgebra added to Project. Check for radial sampling. Need to have Nx, etc. be sequence structure and send separate matrix size for recon parameters. modified: Project.toml modified: src/KomaMRICore.jl modified: src/rawdata/ISMRMRD.jl modified: ../Project.toml --- KomaMRICore/Project.toml | 1 + KomaMRICore/src/KomaMRICore.jl | 1 + KomaMRICore/src/rawdata/ISMRMRD.jl | 22 +++++++++++++--------- Project.toml | 2 ++ 4 files changed, 17 insertions(+), 9 deletions(-) diff --git a/KomaMRICore/Project.toml b/KomaMRICore/Project.toml index 577c46f71..bcf7629c6 100644 --- a/KomaMRICore/Project.toml +++ b/KomaMRICore/Project.toml @@ -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" diff --git a/KomaMRICore/src/KomaMRICore.jl b/KomaMRICore/src/KomaMRICore.jl index f01c8bd25..e8c807659 100644 --- a/KomaMRICore/src/KomaMRICore.jl +++ b/KomaMRICore/src/KomaMRICore.jl @@ -22,6 +22,7 @@ include("simulation/Functors.jl") include("simulation/SimulatorCore.jl") # ISMRMRD +using LinearAlgebra #export signal_to_raw_data export signal_to_raw_data, estimate_seq_recon_dimension # *** CAC 240708 for debugging # Simulator diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 45b4e9ea7..0841a8bb1 100644 --- a/KomaMRICore/src/rawdata/ISMRMRD.jl +++ b/KomaMRICore/src/rawdata/ISMRMRD.jl @@ -282,9 +282,11 @@ function estimate_seq_recon_dimension(seq; sim_params=Dict{String,Any}(), 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, need to & with ADC *** CAC 240708 + 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" @@ -310,7 +312,7 @@ function estimate_seq_recon_dimension(seq; sim_params=Dict{String,Any}(), end if Np_seq > 5*Nv_seq seq_spiral = true - elseif Ns_seq == Nv_seq + elseif sum( k_radius_norm .< 3/Np_seq) >= Nv_seq #Center of K-space highly sampled seq_radial=true else seq_cartesean=true @@ -330,7 +332,9 @@ function estimate_seq_recon_dimension(seq; sim_params=Dict{String,Any}(), Nz = Int64(get(seq.DEF, "Nz", ceil(Int64, Nv_seq/Ny_k))) #pe2 end else - Ny = ceil(Int64, Nv_seq/Ns) + @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 @@ -346,11 +350,11 @@ function estimate_seq_recon_dimension(seq; sim_params=Dict{String,Any}(), if haskey(seq.DEF, "FOV") #needs info or warning for other substitutions??? FOVx, FOVy, FOVz = seq.DEF["FOV"] #[m] - if FOVx > 1.0 - FOVx *= 1e-3 - FOVy *= 1e-3 - FOVz *= 1e-3 - @warn "Scaling FOV to m from older pulseq file mm." + 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 diff --git a/Project.toml b/Project.toml index bfb933857..1a42078f3 100644 --- a/Project.toml +++ b/Project.toml @@ -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"