diff --git a/KomaMRIBase/src/datatypes/Sequence.jl b/KomaMRIBase/src/datatypes/Sequence.jl index febad1c86..ddec4a5d9 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 diff --git a/KomaMRICore/Project.toml b/KomaMRICore/Project.toml index 8be0ccb09..e3d8c42d7 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 3486329ed..6a91080b6 100644 --- a/KomaMRICore/src/KomaMRICore.jl +++ b/KomaMRICore/src/KomaMRICore.jl @@ -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 diff --git a/KomaMRICore/src/rawdata/ISMRMRD.jl b/KomaMRICore/src/rawdata/ISMRMRD.jl index 1ee1ae048..55b70cf75 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=2, use_ndseq=false) Transforms the raw signal into a RawAcquisitionData struct (nearly equivalent to the ISMRMRD format) used for reconstruction with MRIReco. @@ -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 @@ -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[:])) @@ -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 @@ -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] @@ -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( @@ -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 @@ -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) diff --git a/KomaMRICore/src/simulation/GPUFunctions.jl b/KomaMRICore/src/simulation/GPUFunctions.jl index 75e63dbb1..9ef13fe14 100644 --- a/KomaMRICore/src/simulation/GPUFunctions.jl +++ b/KomaMRICore/src/simulation/GPUFunctions.jl @@ -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 @@ -97,4 +96,4 @@ function print_devices(use_gpu = true) _print_devices(backend) end -export print_devices, set_device! +export print_devices 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 diff --git a/Project.toml b/Project.toml index 64e04561f..851c99d1a 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" diff --git a/src/KomaUI.jl b/src/KomaUI.jl index 72ec17e29..a115b4e13 100644 --- a/src/KomaUI.jl +++ b/src/KomaUI.jl @@ -219,19 +219,27 @@ 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 + 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 = raw_aux.params["reconSize"][1:2] - rec_params[:reconSize] = (Nx, Ny) 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 + 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 - @info "Running reconstruction ..." 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!") diff --git a/src/ui/ExportUIFunctions.jl b/src/ui/ExportUIFunctions.jl index 071871943..e2c0746ce 100644 --- a/src/ui/ExportUIFunctions.jl +++ b/src/ui/ExportUIFunctions.jl @@ -123,7 +123,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(