WIP: Implement partial reading for ISMRMRD files

.github/workflows/ci.yml

@@ -45,8 +45,8 @@ jobs:
  shell: julia --color=yes --project=. {0}
  run: |
  using Pkg
- Pkg.add(PackageSpec(path=pwd(), subdir="MRIBase", rev="master"))
- Pkg.add(PackageSpec(path=pwd(), subdir="MRIFiles", rev="master"))
+ Pkg.add(PackageSpec(path=pwd(), subdir="MRIBase"))
+ Pkg.add(PackageSpec(path=pwd(), subdir="MRIFiles"))
  - uses: julia-actions/julia-buildpkg@v1
  - uses: julia-actions/julia-runtest@v1
  - uses: julia-actions/julia-processcoverage@v1
@@ -65,8 +65,8 @@ jobs:
  shell: julia --color=yes --project=. {0}
  run: |
  using Pkg
- Pkg.add(PackageSpec(path=pwd(), subdir="MRIBase", rev="master"))
- Pkg.add(PackageSpec(path=pwd(), subdir="MRIFiles", rev="master"))
+ Pkg.add(PackageSpec(path=pwd(), subdir="MRIBase"))
+ Pkg.add(PackageSpec(path=pwd(), subdir="MRIFiles"))
  - run: |
  julia --project=docs -e '
  using Pkg

MRIFiles/src/ISMRMRD/HDF5Types.jl

@@ -219,6 +219,10 @@ struct HDF5_Acquisition
  data::HVL_T{Float32} #for some reason this is stored as float
 end
 
+struct HDF5_AcquisitionHeaderOnly
+ head::AcquisitionHeaderImmutable
+end
+
 function get_hdf5type_acquisition()
  off = i -> Cint(fieldoffset(HDF5_Acquisition,i))
  datatype = h5t_create(H5T_COMPOUND, sizeof(HDF5_Acquisition) )
@@ -240,6 +244,16 @@ function get_hdf5type_acquisition()
  return datatype
 end
 
+function get_hdf5type_acquisition_header_only()
+ off = i -> Cint(fieldoffset(HDF5_AcquisitionHeaderOnly,i))
+ datatype = h5t_create(H5T_COMPOUND, sizeof(HDF5_AcquisitionHeaderOnly) )
+ d1 = get_hdf5type_acquisitionheader()
+ h5t_insert(datatype, "head", off(1), d1)
+ h5t_close(d1)
+
+ return datatype
+end
+
 
 function writeProfiles(file, dataset, profiles::Vector{Profile})
 

MRIFiles/src/ISMRMRD/ISMRMRD.jl

@@ -12,30 +12,44 @@ end
 
 reads the `ISMRMRDFile` f and stores the result in a `RawAcquisitionDataObject`
 """
-function MRIBase.RawAcquisitionData(f::ISMRMRDFile, dataset="dataset")
+function MRIBase.RawAcquisitionData(f::ISMRMRDFile, dataset="dataset"; 
+ slice=nothing, repetition=nothing, contrast=nothing)
 
  h5open(f.filename) do h
  headerStr = read(h["/$(dataset)/xml"])
  params = GeneralParameters(headerStr[1])
 
- d = read(h["/$(dataset)/data"])
- M = length(d)
+ M = size(h["/$(dataset)/data"]) # for some reason this is a matrix 
+
+ # In the next lines we only read the headers of the profiles such that
+ # we can later filter on them
+ dTypeHeader = get_hdf5type_acquisition_header_only()
+ headerData = Array{HDF5_AcquisitionHeaderOnly}(undef, M)
+ HDF5.API.h5d_read(h["/$(dataset)/data"], dTypeHeader, H5S_ALL, H5S_ALL, H5P_DEFAULT, headerData)
+
  profiles = Profile[]
 
- for m=1:M
- head = read_header(d[m].head)
- D = Int(head.trajectory_dimensions)
- chan = Int(head.active_channels)
+ for m in CartesianIndices(M)
 
- traj = isempty(d[m].traj) ? Matrix{Float32}(undef,0,0) : reshape(d[m].traj, D, :)
+ if (repetition == nothing || headerData[m].head.idx.repetition in repetition) &&
+ (slice == nothing || headerData[m].head.idx.slice in slice) &&
+ (contrast == nothing || headerData[m].head.idx.contrast in contrast)
 
- if !isempty(d[m].data)
- dat = reshape(reinterpret(ComplexF32,d[m].data), :, chan)
- else
- dat = Matrix{ComplexF32}(undef,0,0)
- end
+ d = h["/$(dataset)/data"][m] # Here we read the header again. This could/should be avoided
+ head = read_header(d.head)
+ D = Int(head.trajectory_dimensions)
+ chan = Int(head.active_channels)
 
- push!(profiles, Profile(head,traj,dat) )
+ traj = isempty(d.traj) ? Matrix{Float32}(undef,0,0) : reshape(d.traj, D, :)
+
+ if !isempty(d.data)
+ dat = reshape(reinterpret(ComplexF32,d.data), :, chan)
+ else
+ dat = Matrix{ComplexF32}(undef,0,0)
+ end
+
+ push!(profiles, Profile(head,traj,dat) )
+ end
  end
  return RawAcquisitionData(params, profiles)
  end

test/testISMRMRD.jl

@@ -1,6 +1,5 @@
 using MRIReco
 using Test
-using HTTP
 
 @testset "ISMRMRD" begin
 
@@ -37,6 +36,14 @@ IrecoCopy = reconstruction(AcquisitionData(acqCopy), params)
 
 @test IrecoCopy == Ireco
 
+### test partial read
+
+acqPartial = RawAcquisitionData(f, slice=1) # slice 1 is not contained in the file
+
+@test length(acqPartial.profiles) == 0
+
+### test spiral data
+
 filename = joinpath(datadir, "simple_spiral.h5")
 
 f = ISMRMRDFile(filename)