diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 9a7b759c..d211b3fe 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -17,6 +17,7 @@ jobs: matrix: version: - '1' # Leave this line unchanged. '1' will automatically expand to the latest stable 1.x release of Julia. + - '1.10' #LTS #- 'nightly' os: - ubuntu-latest diff --git a/MRIBase/src/Datatypes/Flags.jl b/MRIBase/src/Datatypes/Flags.jl new file mode 100644 index 00000000..5b39162a --- /dev/null +++ b/MRIBase/src/Datatypes/Flags.jl @@ -0,0 +1,200 @@ +export FLAGS, create_flag_bitmask, flag_is_set, flag_set!, flag_remove!, flag_remove_all!, flags_of + +""" + FLAGS::Dict{String, Int} + +A dictionary mapping string keys (representing flag names) to bitmask values. +Flags are used to indicate specific attributes of the corresponding Profile. Example flag names include "ACQ_FIRST_IN_ENCODE_STEP1", "ACQ_LAST_IN_SLICE", etc. +""" +FLAGS = Dict( + "ACQ_FIRST_IN_ENCODE_STEP1" => 1, + "ACQ_LAST_IN_ENCODE_STEP1" => 2, + "ACQ_FIRST_IN_ENCODE_STEP2" => 3, + "ACQ_LAST_IN_ENCODE_STEP2" => 4, + "ACQ_FIRST_IN_AVERAGE" => 5, + "ACQ_LAST_IN_AVERAGE" => 6, + "ACQ_FIRST_IN_SLICE" => 7, + "ACQ_LAST_IN_SLICE" => 8, + "ACQ_FIRST_IN_CONTRAST" => 9, + "ACQ_LAST_IN_CONTRAST" => 10, + "ACQ_FIRST_IN_PHASE" => 11, + "ACQ_LAST_IN_PHASE" => 12, + "ACQ_FIRST_IN_REPETITION" => 13, + "ACQ_LAST_IN_REPETITION" => 14, + "ACQ_FIRST_IN_SET" => 15, + "ACQ_LAST_IN_SET" => 16, + "ACQ_FIRST_IN_SEGMENT" => 17, + "ACQ_LAST_IN_SEGMENT" => 18, + "ACQ_IS_NOISE_MEASUREMENT" => 19, + "ACQ_IS_PARALLEL_CALIBRATION" => 20, + "ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING" => 21, + "ACQ_IS_REVERSE" => 22, + "ACQ_IS_NAVIGATION_DATA" => 23, + "ACQ_IS_PHASECORR_DATA" => 24, + "ACQ_LAST_IN_MEASUREMENT" => 25, + "ACQ_IS_HPFEEDBACK_DATA" => 26, + "ACQ_IS_DUMMYSCAN_DATA" => 27, + "ACQ_IS_RTFEEDBACK_DATA" => 28, + "ACQ_IS_SURFACECOILCORRECTIONSCAN_DATA" => 29, + "ACQ_COMPRESSION1" => 53, + "ACQ_COMPRESSION2" => 54, + "ACQ_COMPRESSION3" => 55, + "ACQ_COMPRESSION4" => 56, + "ACQ_USER1" => 57, + "ACQ_USER2" => 58, + "ACQ_USER3" => 59, + "ACQ_USER4" => 60, + "ACQ_USER5" => 61, + "ACQ_USER6" => 62, + "ACQ_USER7" => 63, + "ACQ_USER8" => 64 +) + +function bitshift(A, k) + k >= 0 ? out = A << k : out = A >> abs(k) + return out +end + +""" + create_flag_bitmask(flag::AbstractString) -> UInt64 + create_flag_bitmask(flag::Integer) -> UInt64 + +Creates a bitmask for the given flag. + +- If the flag is a string, its integer value is looked up in `FLAGS`. +- If the flag is an integer, it must be positive. The bitmask is created by left-shifting `1` to the corresponding bit position. + +# Arguments +- `flag`: A string representing the flag name or a positive integer. + +# Returns +- A `UInt64` bitmask where only the bit corresponding to the flag is set. + +# Throws +- `KeyError` if the string flag is not found in `FLAGS`. +- `DomainError` if the integer flag is not positive. + +# Example +```julia +create_flag_bitmask("ACQ_FIRST_IN_ENCODE_STEP1") # 0x0000000000000001 +create_flag_bitmask(2) # 0x0000000000000002 +``` +""" +create_flag_bitmask(flag::T) where T = error("Unexpected type for bitmask, expected String or positive Integer, found $T") +create_flag_bitmask(flag::AbstractString) = create_flag_bitmask(FLAGS[flag]) +function create_flag_bitmask(flag::Integer) + flag > 0 || throw(DomainError(flag, "Bitmask can only be created for positive integers")) + b = UInt64(flag) + return bitshift(UInt64(1), b - 1) +end + +""" + flag_is_set(obj::Profile, flag) -> Bool + +Checks whether a specific flag is set in the given `Profile` object. + +# Arguments +- `obj`: A `Profile` object with a `head` field containing a `flags` attribute (as a bitmask). +- `flag`: A string (flag name) or an integer (bit position). + +# Returns +- `true` if the flag is set in the `obj.head.flags` bitmask; otherwise, `false`. + +# Example +```julia +flag_is_set(profile, "ACQ_FIRST_IN_ENCODE_STEP1") # true or false +``` +""" +function flag_is_set(obj::Profile, flag) + bitmask = create_flag_bitmask(flag) + ret = obj.head.flags & bitmask > 0 + return ret +end + +""" + flag_set!(obj::Profile, flag) + +Sets the specified flag in the given `Profile` object. + +# Arguments +- `obj`: A `Profile` object with a `head` field containing a `flags` attribute (as a bitmask). +- `flag`: A string (flag name) or an integer (bit position). + +# Modifies +- `obj.head.flags` by setting the bit corresponding to the flag. + +# Example +```julia +flag_set!(profile, "ACQ_FIRST_IN_SLICE") +``` +""" +function flag_set!(obj::Profile, flag) + bitmask = create_flag_bitmask(flag) + obj.head.flags = obj.head.flags | bitmask +end + +""" + flag_remove!(obj::Profile, flag) + +Removes (clears) the specified flag in the given `Profile` object. + +# Arguments +- `obj`: A `Profile` object with a `head` field containing a `flags` attribute (as a bitmask). +- `flag`: A string (flag name) or an integer (bit position). + +# Modifies +- `obj.head.flags` by clearing the bit corresponding to the flag. + +# Example +```julia +flag_remove!(profile, "ACQ_LAST_IN_PHASE") +``` +""" +function flag_remove!(obj::Profile, flag) + bitmask = create_flag_bitmask(flag) + obj.head.flags = obj.head.flags & ~bitmask +end + +""" + flag_remove_all!(obj::Profile) + +Clears all flags in the given `Profile` object. + +# Arguments +- `obj`: A `Profile` object with a `head` field containing a `flags` attribute (as a bitmask). + +# Modifies +- `obj.head.flags`, setting it to `0` (all flags cleared). + +# Example +```julia +flag_remove_all!(profile) +``` +""" +function flag_remove_all!(obj::Profile) + obj.head.flags = UInt64(0); +end + +""" + flags_of(obj::Profile) -> Vector{String} + +Returns a list of all flags that are set in the given `Profile` object. + +# Arguments +- `obj`: A `Profile` object with a `head` field containing a `flags` attribute (as a bitmask). + +# Returns +- An array of strings representing the names of the flags that are currently set. + +# Example +```julia +flags = flags_of(profile) # ["ACQ_FIRST_IN_ENCODE_STEP1", "ACQ_LAST_IN_SLICE"] +``` +""" +function flags_of(obj::Profile) + flags_ = String[] + for f in FLAGS + flag_is_set(obj::Profile, f.first) ? push!(flags_,f.first) : nothing + end + return flags_ +end \ No newline at end of file diff --git a/MRIBase/src/Datatypes/RawAcqData.jl b/MRIBase/src/Datatypes/RawAcqData.jl index d0491d29..ac4dee17 100644 --- a/MRIBase/src/Datatypes/RawAcqData.jl +++ b/MRIBase/src/Datatypes/RawAcqData.jl @@ -206,7 +206,7 @@ function subsampleIndices(f::RawAcquisitionData; slice::Int=1, contrast::Int=1, i1 = tr_center_idx - center_sample + 1 + f.profiles[i].head.discard_pre i2 = tr_center_idx - center_sample + numSamp - f.profiles[i].head.discard_post # convert to linear index - lineIdx = collect(i1:i2) .+ numSamp*((encSt2[i]-1)*numProf + (encSt1[i]-1)) + lineIdx = collect(i1:i2) .+ numEncSamp *((encSt2[i]-1)*numProf + (encSt1[i]-1)) append!(idx, lineIdx) end @@ -264,8 +264,9 @@ function rawdata(f::RawAcquisitionData; slice::Int=1, contrast::Int=1, repetitio i1 = f.profiles[l].head.discard_pre + 1 i2 = i1+numSampPerProfile-1 kdata[:,cnt,:] .= f.profiles[l].data[i1:i2, :] - if f.profiles[l].head.read_dir[1] < 0 + if flag_is_set(f.profiles[l],"ACQ_IS_REVERSE") kdata[:,cnt,:] .= reverse(kdata[:,cnt,:], dims=1) + flag_remove!(f.profiles[l],"ACQ_IS_REVERSE") end cnt += 1 end diff --git a/MRIBase/src/MRIBase.jl b/MRIBase/src/MRIBase.jl index 9b262137..73348c06 100644 --- a/MRIBase/src/MRIBase.jl +++ b/MRIBase/src/MRIBase.jl @@ -5,6 +5,5 @@ using NFFTTools # for density compensation weights in trajectory include("Trajectories/Trajectories.jl") include("Datatypes/Datatypes.jl") - - +include("Datatypes/Flags.jl") end # module diff --git a/MRIBase/test/runtests.jl b/MRIBase/test/runtests.jl index 18d41534..e2cca7b0 100644 --- a/MRIBase/test/runtests.jl +++ b/MRIBase/test/runtests.jl @@ -1,3 +1,4 @@ using Test, MRIBase -include("testTrajectories.jl") \ No newline at end of file +include("testTrajectories.jl") +include("testFlags.jl") \ No newline at end of file diff --git a/MRIBase/test/testFlags.jl b/MRIBase/test/testFlags.jl new file mode 100644 index 00000000..2dd971fa --- /dev/null +++ b/MRIBase/test/testFlags.jl @@ -0,0 +1,34 @@ + +@testset "Flags" begin + head = AcquisitionHeader() + traj = Matrix{Float32}(undef,0,0) + dat = Matrix{ComplexF32}(undef,0,0) + p = Profile(head,traj,dat) + + p2 = deepcopy(p) + @test ~flag_is_set(p,"ACQ_IS_REVERSE") + flag_set!(p,"ACQ_IS_REVERSE") + flag_set!(p2,22) + @test flag_is_set(p,"ACQ_IS_REVERSE") + @test p2.head.flags == p.head.flags + flag_remove!(p,"ACQ_IS_REVERSE") + @test ~flag_is_set(p,22) + + flag_set!(p,"ACQ_IS_REVERSE") + flag_set!(p,"ACQ_IS_NAVIGATION_DATA") + + fl = flags_of(p) + + @test isequal(fl[1],"ACQ_IS_NAVIGATION_DATA") && isequal(fl[2],"ACQ_IS_REVERSE") + @test flag_is_set(p,"ACQ_IS_REVERSE") & flag_is_set(p,"ACQ_IS_NAVIGATION_DATA") + flag_remove_all!(p) + @test p.head.flags == UInt64(0) + + @test_throws Exception flag_set!(p, -1) + @test_throws Exception flag_set!(p, "NOT_A_VALID_FLAG") + @test_throws Exception flag_is_set(p, -1) + @test_throws Exception flag_is_set(p, "NOT_A_VALID_FLAG") + @test_throws Exception flag_remove!(p, -1) + @test_throws Exception flag_remove!(p, "NOT_A_VALID_FLAG") + end + \ No newline at end of file diff --git a/MRIFiles/src/Bruker/Bruker.jl b/MRIFiles/src/Bruker/Bruker.jl index 503e2309..25653010 100644 --- a/MRIFiles/src/Bruker/Bruker.jl +++ b/MRIFiles/src/Bruker/Bruker.jl @@ -29,7 +29,7 @@ mutable struct BrukerFile <: MRIFile maxEntriesAcqp end -include("BrukerSequence.jl") +include("BrukerSequencePV6.jl") include("BrukerSequence360.jl") function BrukerFile(path::String; maxEntriesAcqp=2000) @@ -269,135 +269,6 @@ function MRIBase.RawAcquisitionData(b::BrukerFile) end -function RawAcquisitionDataRawDataSpiral(b::BrukerFile) - # Have not a way to read out this. Not sure if its true - T = Complex{Int32} - - filename = joinpath(b.path, "rawdata.job0") - - profileLength = pvmTrajSamples(b) - phaseFactor = acqPhaseFactor(b) - numSlices = acqNumSlices(b) - numEchos = acqNumEchos(b) - numRep = acqNumRepetitions(b) - numProfiles = pvmTrajInterleaves(b) - numChannels = acqNumCoils(b) - - I = open(filename,"r") do fd - read!(fd,Array{T,7}(undef, profileLength, - numChannels, - numEchos, - phaseFactor, - numSlices, - numProfiles, - numRep)) - end - - # TODO the following of course is just for debugging - return I -end - -function RawAcquisitionDataFid(b::BrukerFile) - T = Complex{acqDataType(b)} - - filename = joinpath(b.path, "fid") - - N = acqSize(b) - # The data is padded in case it is not a multiple of 1024 - # For multi-channel acquisition data at concatenate then padded to a multiple of 1024 bytes - numChannel = pvmEncNReceivers(b) - profileLength = Int((ceil(N[1]*numChannel*sizeof(T)/1024))*1024/sizeof(T)) - numAvailableChannel = pvmEncAvailReceivers(b) - phaseFactor = acqPhaseFactor(b) - numSlices = acqNumSlices(b) - numEchos = acqNumEchos(b) - numEncSteps2 = length(N) == 3 ? N[3] : 1 - numRep = acqNumRepetitions(b) - - I = open(filename,"r") do fd - read!(fd,Array{T,7}(undef, profileLength, - numEchos, - phaseFactor, - numSlices, - div(N[2], phaseFactor), - numEncSteps2, - numRep))[1:N[1]*numChannel,:,:,:,:,:,:] - end - - encSteps1 = pvmEncSteps1(b) - encSteps1 = encSteps1.-minimum(encSteps1) - - encSteps2 = pvmEncSteps2(b) - encSteps2 = encSteps2.-minimum(encSteps2) - - objOrd = acqObjOrder(b) - objOrd = objOrd.-minimum(objOrd) - - gradMatrix = acqGradMatrix(b) - - offset1 = acqReadOffset(b) - offset2 = acqPhase1Offset(b) - offset3 = ndims(b) == 2 ? acqSliceOffset(b) : pvmEffPhase2Offset(b) - - profiles = Profile[] - for nR = 1:numRep - for nEnc2 = 1:numEncSteps2 - for nPhase2 = 1:div(N[2], phaseFactor) - for nSl = 1:numSlices - for nPhase1 = 1:phaseFactor - for nEcho=1:numEchos - counter = EncodingCounters(kspace_encode_step_1=encSteps1[nPhase1+phaseFactor*(nPhase2-1)], - kspace_encode_step_2=encSteps2[nEnc2], - average=0, - slice=objOrd[nSl], - contrast=nEcho-1, - phase=0, - repetition=nR-1, - set=0, - segment=0 ) - - G = gradMatrix[:,:,nSl] - read_dir = (G[1,1],G[2,1],G[3,1]) - phase_dir = (G[1,2],G[2,2],G[3,2]) - slice_dir = (G[1,3],G[2,3],G[3,3]) - - # Not sure if the following is correct... - pos = offset1[nSl]*G[:,1] + - offset2[nSl]*G[:,2] + - offset3[nSl]*G[:,3] - - position = (pos[1], pos[2], pos[3]) - - head = AcquisitionHeader(number_of_samples=N[1], idx=counter, - read_dir=read_dir.*(-1)^(nEcho-1), phase_dir=phase_dir, - slice_dir=slice_dir, position=position, - center_sample=div(N[1],2), - available_channels = numAvailableChannel, #TODO - active_channels = numChannel) - traj = Matrix{Float32}(undef,0,0) - dat = map(T, reshape(I[:,nEcho,nPhase1,nSl,nPhase2,nEnc2,nR],:,numChannel)) - push!(profiles, Profile(head,traj,dat) ) - end - end - end - end - end - end - - params = brukerParams(b) - params["trajectory"] = "cartesian" - N = acqSize(b) - if length(N) < 3 - N_ = ones(Int,3) - N_[1:length(N)] .= N - N = N_ - end - params["encodedSize"] = N - - return RawAcquisitionData(params, profiles) -end - - ##### Reco function recoData(f::BrukerFile,procno::Int=1) recoFilename = joinpath(f.path,"pdata", string(procno), "2dseq") diff --git a/MRIFiles/src/Bruker/BrukerSequence.jl b/MRIFiles/src/Bruker/BrukerSequence.jl deleted file mode 100644 index ca32f027..00000000 --- a/MRIFiles/src/Bruker/BrukerSequence.jl +++ /dev/null @@ -1,97 +0,0 @@ -export RawAcquisitionData_3DUTE - -function RawAcquisitionData_3DUTE(b::BrukerFile) - T = Complex{acqDataType(b)} - - filename = joinpath(b.path, "fid") - filenameTraj = joinpath(b.path, "traj") - - N = acqSize(b) - # The data is padded in case it is not a multiple of 1024 - # For multi-channel acquisition data at concatenate then padded to a multiple of 1024 bytes - numChannel = pvmEncNReceivers(b) - profileLength = Int((ceil(N[1]*numChannel*sizeof(T)/1024))*1024/sizeof(T)) - numAvailableChannel = pvmEncAvailReceivers(b) - phaseFactor = acqPhaseFactor(b) - numSlices = acqNumSlices(b) - numEchos = acqNumEchos(b) - numEncSteps2 = length(N) == 3 ? N[3] : 1 - numRep = acqNumRepetitions(b) - sample_time = parse(Float64,b["PVM_DigDw"])/1000 #us - - I = open(filename,"r") do fd - read!(fd,Array{T,7}(undef, profileLength, - numEchos, - phaseFactor, - numSlices, - div(N[2], phaseFactor), - numEncSteps2, - numRep))[1:N[1]*numChannel,:,:,:,:,:,:] - end - - traj = open(filenameTraj,"r") do fd - read!(fd,Array{Float64,3}(undef, 3,N[1], - N[2])) - end - - objOrd = acqObjOrder(b) - objOrd = objOrd.-minimum(objOrd) - - gradMatrix = acqGradMatrix(b) - - offset1 = acqReadOffset(b) - offset2 = acqPhase1Offset(b) - offset3 = ndims(b) == 2 ? acqSliceOffset(b) : pvmEffPhase2Offset(b) - - profiles = Profile[] - for nR = 1:numRep - for nEnc2 = 1:numEncSteps2 - for nPhase2 = 1:div(N[2], phaseFactor) - for nSl = 1:numSlices - for nPhase1 = 1:phaseFactor - for nEcho=1:numEchos - counter = EncodingCounters(kspace_encode_step_1=0, - kspace_encode_step_2=0, - average=0, - slice=objOrd[nSl], - contrast=nEcho-1, - phase=0, - repetition=nR-1, - set=0, - segment=0 ) - - G = gradMatrix[:,:,nSl] - read_dir = (G[1,1],G[2,1],G[3,1]) - phase_dir = (G[1,2],G[2,2],G[3,2]) - slice_dir = (G[1,3],G[2,3],G[3,3]) - - # Not sure if the following is correct... - pos = offset1[nSl]*G[:,1] + - offset2[nSl]*G[:,2] + - offset3[nSl]*G[:,3] - - position = (pos[1], pos[2], pos[3]) - - head = AcquisitionHeader(number_of_samples=N[1], idx=counter, - read_dir=read_dir, phase_dir=phase_dir, - slice_dir=slice_dir, position=position, - center_sample=0, - available_channels = numAvailableChannel, #TODO - active_channels = numChannel, - trajectory_dimensions = 3, - sample_time_us = sample_time) - - dat = map(T, reshape(I[:,nEcho,nPhase1,nSl,nPhase2,nEnc2,nR],:,numChannel)) - push!(profiles, Profile(head,convert.(Float32,traj[:,:,nPhase2]),dat) ) - end - end - end - end - end - end - - params = brukerParams(b) - params["trajectory"] = "custom" - params["encodedSize"] = pvmMatrix(b) - return RawAcquisitionData(params, profiles) -end diff --git a/MRIFiles/src/Bruker/BrukerSequence360.jl b/MRIFiles/src/Bruker/BrukerSequence360.jl index d876a1b0..6a39a8e0 100644 --- a/MRIFiles/src/Bruker/BrukerSequence360.jl +++ b/MRIFiles/src/Bruker/BrukerSequence360.jl @@ -21,9 +21,11 @@ function RawAcquisitionDataFid_360(b::BrukerFile) numEchos = MRIFiles.acqNumEchos(b) phaseFactor = MRIFiles.acqPhaseFactor(b) numRep = MRIFiles.acqNumRepetitions(b) - numEncSteps = MRIFiles.acqSpatialSize1(b) + numEncSteps = MRIFiles.acqSpatialSize1(b) + readoutLength = parse.(Int,b["PVM_EncMatrix"])[1] - profileLength = N[1]*numChannel#MRIFiles.acqSize(b)[1] #Int((ceil(N[1]*numChannel*sizeof(dtype)/1024))*1024/sizeof(dtype)) # number of points + zeros + centerSample = parse.(Int,b["PVM_EncPftOverscans"])[1] + profileLength = readoutLength*numChannel#MRIFiles.acqSize(b)[1] #Int((ceil(N[1]*numChannel*sizeof(dtype)/1024))*1024/sizeof(dtype)) # number of points + zeros I = open(filename,"r") do fd read!(fd,Array{T,6}(undef, profileLength, @@ -34,8 +36,8 @@ function RawAcquisitionDataFid_360(b::BrukerFile) numRep)) end - encSteps1 = parse.(Int,b["PVM_EncGenSteps1"]).+round(Int,N[2]/2) - encSteps2 = parse.(Int,b["PVM_EncGenSteps2"]).+round(Int,N[3]/2) + encSteps1 = parse.(Int,b["PVM_EncGenSteps1"]).+floor(Int,N[2]/2) + encSteps2 = parse.(Int,b["PVM_EncGenSteps2"]).+floor(Int,N[3]/2) objOrd = MRIFiles.acqObjOrder(b) objOrd = objOrd.-minimum(objOrd) @@ -47,7 +49,7 @@ function RawAcquisitionDataFid_360(b::BrukerFile) offset3 = ndims(b) == 2 ? MRIFiles.acqSliceOffset(b) : MRIFiles.pvmEffPhase2Offset(b) - + profiles = Profile[] for nR = 1:numRep for nEnc = 1:div(numEncSteps, phaseFactor) @@ -79,12 +81,23 @@ function RawAcquisitionDataFid_360(b::BrukerFile) head = AcquisitionHeader(number_of_samples=N[1], idx=counter, read_dir=read_dir, phase_dir=phase_dir, slice_dir=slice_dir, position=position, - center_sample=div(N[1],2), + center_sample=centerSample,#div(N[1],2), available_channels = numChannel, #numAvailableChannel ? active_channels = numChannel) + + traj = Matrix{Float32}(undef,0,0) dat = map(T, reshape(I[:,nEcho,nPhase,nSl,nEnc,nR],:,numChannel)) - push!(profiles, Profile(head,traj,dat) ) + p = Profile(head,traj,dat) + + if (b["EchoAcqMode"] == "allEchoes") && mod(nEcho,2) == 0 # set reverse specific flags for MGE + @info nEcho + flag_set!(p,"ACQ_IS_REVERSE") + end + + push!(profiles,p) + + end end end diff --git a/MRIFiles/src/Bruker/BrukerSequencePV6.jl b/MRIFiles/src/Bruker/BrukerSequencePV6.jl new file mode 100644 index 00000000..c6089099 --- /dev/null +++ b/MRIFiles/src/Bruker/BrukerSequencePV6.jl @@ -0,0 +1,234 @@ +export RawAcquisitionData_3DUTE + + +function RawAcquisitionDataRawDataSpiral(b::BrukerFile) + # Have not a way to read out this. Not sure if its true + T = Complex{Int32} + + filename = joinpath(b.path, "rawdata.job0") + + profileLength = pvmTrajSamples(b) + phaseFactor = acqPhaseFactor(b) + numSlices = acqNumSlices(b) + numEchos = acqNumEchos(b) + numRep = acqNumRepetitions(b) + numProfiles = pvmTrajInterleaves(b) + numChannels = acqNumCoils(b) + + I = open(filename,"r") do fd + read!(fd,Array{T,7}(undef, profileLength, + numChannels, + numEchos, + phaseFactor, + numSlices, + numProfiles, + numRep)) + end + + # TODO the following of course is just for debugging + return I +end + +function RawAcquisitionDataFid(b::BrukerFile) + T = Complex{acqDataType(b)} + + filename = joinpath(b.path, "fid") + + N = acqSize(b) + # The data is padded in case it is not a multiple of 1024 + # For multi-channel acquisition data at concatenate then padded to a multiple of 1024 bytes + numChannel = pvmEncNReceivers(b) + profileLength = Int((ceil(N[1]*numChannel*sizeof(T)/1024))*1024/sizeof(T)) + numAvailableChannel = pvmEncAvailReceivers(b) + phaseFactor = acqPhaseFactor(b) + numSlices = acqNumSlices(b) + numEchos = acqNumEchos(b) + numEncSteps2 = length(N) == 3 ? N[3] : 1 + numRep = acqNumRepetitions(b) + + I = open(filename,"r") do fd + read!(fd,Array{T,7}(undef, profileLength, + numEchos, + phaseFactor, + numSlices, + div(N[2], phaseFactor), + numEncSteps2, + numRep))[1:N[1]*numChannel,:,:,:,:,:,:] + end + + encSteps1 = pvmEncSteps1(b) + encSteps1 = encSteps1.-minimum(encSteps1) + + encSteps2 = pvmEncSteps2(b) + encSteps2 = encSteps2.-minimum(encSteps2) + + objOrd = acqObjOrder(b) + objOrd = objOrd.-minimum(objOrd) + + gradMatrix = acqGradMatrix(b) + + offset1 = acqReadOffset(b) + offset2 = acqPhase1Offset(b) + offset3 = ndims(b) == 2 ? acqSliceOffset(b) : pvmEffPhase2Offset(b) + + profiles = Profile[] + for nR = 1:numRep + for nEnc2 = 1:numEncSteps2 + for nPhase2 = 1:div(N[2], phaseFactor) + for nSl = 1:numSlices + for nPhase1 = 1:phaseFactor + for nEcho=1:numEchos + counter = EncodingCounters(kspace_encode_step_1=encSteps1[nPhase1+phaseFactor*(nPhase2-1)], + kspace_encode_step_2=encSteps2[nEnc2], + average=0, + slice=objOrd[nSl], + contrast=nEcho-1, + phase=0, + repetition=nR-1, + set=0, + segment=0 ) + + G = gradMatrix[:,:,nSl] + read_dir = (G[1,1],G[2,1],G[3,1]) + phase_dir = (G[1,2],G[2,2],G[3,2]) + slice_dir = (G[1,3],G[2,3],G[3,3]) + + # Not sure if the following is correct... + pos = offset1[nSl]*G[:,1] + + offset2[nSl]*G[:,2] + + offset3[nSl]*G[:,3] + + position = (pos[1], pos[2], pos[3]) + + head = AcquisitionHeader(number_of_samples=N[1], idx=counter, + read_dir=read_dir.*(-1)^(nEcho-1), phase_dir=phase_dir, + slice_dir=slice_dir, position=position, + center_sample=div(N[1],2), + available_channels = numAvailableChannel, #TODO + active_channels = numChannel) + traj = Matrix{Float32}(undef,0,0) + dat = map(T, reshape(I[:,nEcho,nPhase1,nSl,nPhase2,nEnc2,nR],:,numChannel)) + + p = Profile(head,traj,dat) + if (b["EchoAcqMode"] == "allEchoes") && mod(nEcho,2) == 0 # set reverse specific flags for MGE + @info nEcho + flag_set!(p,"ACQ_IS_REVERSE") + end + + push!(profiles, p ) + end + end + end + end + end + end + + params = brukerParams(b) + params["trajectory"] = "cartesian" + N = acqSize(b) + if length(N) < 3 + N_ = ones(Int,3) + N_[1:length(N)] .= N + N = N_ + end + params["encodedSize"] = N + + return RawAcquisitionData(params, profiles) +end + + +function RawAcquisitionData_3DUTE(b::BrukerFile) + T = Complex{acqDataType(b)} + + filename = joinpath(b.path, "fid") + filenameTraj = joinpath(b.path, "traj") + + N = acqSize(b) + # The data is padded in case it is not a multiple of 1024 + # For multi-channel acquisition data at concatenate then padded to a multiple of 1024 bytes + numChannel = pvmEncNReceivers(b) + profileLength = Int((ceil(N[1]*numChannel*sizeof(T)/1024))*1024/sizeof(T)) + numAvailableChannel = pvmEncAvailReceivers(b) + phaseFactor = acqPhaseFactor(b) + numSlices = acqNumSlices(b) + numEchos = acqNumEchos(b) + numEncSteps2 = length(N) == 3 ? N[3] : 1 + numRep = acqNumRepetitions(b) + sample_time = parse(Float64,b["PVM_DigDw"])/1000 #us + + I = open(filename,"r") do fd + read!(fd,Array{T,7}(undef, profileLength, + numEchos, + phaseFactor, + numSlices, + div(N[2], phaseFactor), + numEncSteps2, + numRep))[1:N[1]*numChannel,:,:,:,:,:,:] + end + + traj = open(filenameTraj,"r") do fd + read!(fd,Array{Float64,3}(undef, 3,N[1], + N[2])) + end + + objOrd = acqObjOrder(b) + objOrd = objOrd.-minimum(objOrd) + + gradMatrix = acqGradMatrix(b) + + offset1 = acqReadOffset(b) + offset2 = acqPhase1Offset(b) + offset3 = ndims(b) == 2 ? acqSliceOffset(b) : pvmEffPhase2Offset(b) + + profiles = Profile[] + for nR = 1:numRep + for nEnc2 = 1:numEncSteps2 + for nPhase2 = 1:div(N[2], phaseFactor) + for nSl = 1:numSlices + for nPhase1 = 1:phaseFactor + for nEcho=1:numEchos + counter = EncodingCounters(kspace_encode_step_1=0, + kspace_encode_step_2=0, + average=0, + slice=objOrd[nSl], + contrast=nEcho-1, + phase=0, + repetition=nR-1, + set=0, + segment=0 ) + + G = gradMatrix[:,:,nSl] + read_dir = (G[1,1],G[2,1],G[3,1]) + phase_dir = (G[1,2],G[2,2],G[3,2]) + slice_dir = (G[1,3],G[2,3],G[3,3]) + + # Not sure if the following is correct... + pos = offset1[nSl]*G[:,1] + + offset2[nSl]*G[:,2] + + offset3[nSl]*G[:,3] + + position = (pos[1], pos[2], pos[3]) + + head = AcquisitionHeader(number_of_samples=N[1], idx=counter, + read_dir=read_dir, phase_dir=phase_dir, + slice_dir=slice_dir, position=position, + center_sample=0, + available_channels = numAvailableChannel, #TODO + active_channels = numChannel, + trajectory_dimensions = 3, + sample_time_us = sample_time) + + dat = map(T, reshape(I[:,nEcho,nPhase1,nSl,nPhase2,nEnc2,nR],:,numChannel)) + push!(profiles, Profile(head,convert.(Float32,traj[:,:,nPhase2]),dat) ) + end + end + end + end + end + end + + params = brukerParams(b) + params["trajectory"] = "custom" + params["encodedSize"] = pvmMatrix(b) + return RawAcquisitionData(params, profiles) +end diff --git a/MRIFiles/src/Bruker/Jcampdx.jl b/MRIFiles/src/Bruker/Jcampdx.jl index 266f9bb1..4d665c0a 100644 --- a/MRIFiles/src/Bruker/Jcampdx.jl +++ b/MRIFiles/src/Bruker/Jcampdx.jl @@ -77,7 +77,7 @@ function read(file::JcampdxFile, stream::IO, keylist::Vector=String[]; maxEntrie if val[1] != '(' file.dict[key] = val else - if val[2] != ' ' + if val[2] != ' ' || val[3] == '<' file.dict[key] = val if val[end] == ')' file.dict[key] = val diff --git a/MRIFiles/test/testISMRMRD.jl b/MRIFiles/test/testISMRMRD.jl index c49ff982..9e43d7bd 100644 --- a/MRIFiles/test/testISMRMRD.jl +++ b/MRIFiles/test/testISMRMRD.jl @@ -33,7 +33,9 @@ io = IOBuffer() write(io, acq.profiles[1].head) ioCopy = IOBuffer() write(ioCopy, acqCopy.profiles[1].head) -@test io.data == ioCopy.data +seekstart(io) +seekstart(ioCopy) +@test read(io) == read(ioCopy) @test acqCopy.profiles[1].data == acq.profiles[1].data # test that the offsets in the header are not padded diff --git a/docs/src/acquisitionData.md b/docs/src/acquisitionData.md index 6daa6386..2618f085 100644 --- a/docs/src/acquisitionData.md +++ b/docs/src/acquisitionData.md @@ -31,6 +31,71 @@ measured after a single excitation during an MRI experiment. It has members `head`, `traj`, and `data`, which exactly resemble the structures specified by the ISMRMRD file format. +AcquisitionHeader has exactly the same structure as ISMRMRD. You can find more information about it [here](https://ismrmrd.readthedocs.io/en/latest/mrd_raw_data.html#acquisitionheader) + +Two fields are especially important in it : +- flags +- idx + +### flags + +The flags field in the AcquisitionHeader is a 64 bit mask that can be used to indicate specific attributes of the corresponding readout. One usage of these flags is to reverse the signal during conversion from RawAcquisitionData to AcquisitionData if the flag "ACQ_IS_REVERSE" is set. + +```julia +FLAGS = Dict( + "ACQ_FIRST_IN_ENCODE_STEP1" => 1, + "ACQ_LAST_IN_ENCODE_STEP1" => 2, + "ACQ_FIRST_IN_ENCODE_STEP2" => 3, + "ACQ_LAST_IN_ENCODE_STEP2" => 4, + "ACQ_FIRST_IN_AVERAGE" => 5, + "ACQ_LAST_IN_AVERAGE" => 6, + "ACQ_FIRST_IN_SLICE" => 7, + "ACQ_LAST_IN_SLICE" => 8, + "ACQ_FIRST_IN_CONTRAST" => 9, + "ACQ_LAST_IN_CONTRAST" => 10, + "ACQ_FIRST_IN_PHASE" => 11, + "ACQ_LAST_IN_PHASE" => 12, + "ACQ_FIRST_IN_REPETITION" => 13, + "ACQ_LAST_IN_REPETITION" => 14, + "ACQ_FIRST_IN_SET" => 15, + "ACQ_LAST_IN_SET" => 16, + "ACQ_FIRST_IN_SEGMENT" => 17, + "ACQ_LAST_IN_SEGMENT" => 18, + "ACQ_IS_NOISE_MEASUREMENT" => 19, + "ACQ_IS_PARALLEL_CALIBRATION" => 20, + "ACQ_IS_PARALLEL_CALIBRATION_AND_IMAGING" => 21, + "ACQ_IS_REVERSE" => 22, + "ACQ_IS_NAVIGATION_DATA" => 23, + "ACQ_IS_PHASECORR_DATA" => 24, + "ACQ_LAST_IN_MEASUREMENT" => 25, + "ACQ_IS_HPFEEDBACK_DATA" => 26, + "ACQ_IS_DUMMYSCAN_DATA" => 27, + "ACQ_IS_RTFEEDBACK_DATA" => 28, + "ACQ_IS_SURFACECOILCORRECTIONSCAN_DATA" => 29, + "ACQ_COMPRESSION1" => 53, + "ACQ_COMPRESSION2" => 54, + "ACQ_COMPRESSION3" => 55, + "ACQ_COMPRESSION4" => 56, + "ACQ_USER1" => 57, + "ACQ_USER2" => 58, + "ACQ_USER3" => 59, + "ACQ_USER4" => 60, + "ACQ_USER5" => 61, + "ACQ_USER6" => 62, + "ACQ_USER7" => 63, + "ACQ_USER8" => 64 +) +``` + +You can check the flags of a profile with `flags_of(p:Profile)` or `flag_is_set` and manipulate them with thus functions : +- `flag_set!(obj::Profile, flag)` +- `flag_remove!(obj::Profile, flag)` +- `flag_remove_all!(obj::Profile)` + +### idx + +MR acquisitions often loop through a set of counters (e.g. phase encodes) in a complete experiment. The following encoding counters are referred to by the idx field in the AcquisitionHeader ([See the ISMRMRD documentation](https://ismrmrd.readthedocs.io/en/latest/mrd_raw_data.html#mrd-encodingcounters)) + ## Preprocessed Data The `RawAcquisitionData` can be preprocessed into a form, which makes it more diff --git a/test/testReconstruction.jl b/test/testReconstruction.jl index 43f94fe6..9cb4d975 100644 --- a/test/testReconstruction.jl +++ b/test/testReconstruction.jl @@ -614,7 +614,7 @@ function testCSReco3d(N=128;arrayType = Array) Ireco = reconstruction(acqData, params) - @test (norm(vec(I)-vec(Ireco))/norm(vec(I))) < 1e-1 + @test (norm(vec(I)-vec(Ireco))/norm(vec(I))) < 2e-1 end function testCSSenseReco3d(N=128;arrayType = Array)