MagneticResonanceImaging / MRIReco.jl

Julia Package for MRI Reconstruction
https://magneticresonanceimaging.github.io/MRIReco.jl/latest/
Other
85 stars 22 forks source link

Partial (or lazy) Loading of RawAcquisitionData #76

Open tknopp opened 2 years ago

tknopp commented 2 years ago

I am splitting of this feature request from #73 since this would be easier to implement for ISMRMRD.

@mrikasper @alexjaffray: Could you come up with a small simple example of what you have in mind here? Ideally you would design a small minimal example that we need for the test suite anyway. The example can be small and simulation based but it should represent your use case well.

From my understanding one wants something like

   RawAcquisitionData(f::ISMRMRDFile, slices, dataset="dataset")

If that is correct we would need to somehow filter profiles in this for loop

https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/master/MRIFiles/src/ISMRMRD/ISMRMRD.jl#L25

All this should IMHO happen while loading the RawDataAcquisition object.

mrikasper commented 2 years ago

We might be able to create a minimum example from this code

https://github.com/mrikasper/julia-recon-advances-in-spiral-fmri

related to the spiral fMRI paper and data we published.

As you see in the main.jl L65, I already tried a crude version of extracting a single slice only and adjusting the header for that. RawDataAcquisition indeed seems to be the right place. It would also allow via writing back to ISMRMRD files to create reduced file versions for later faster loading, wouldn't it?

What should the minimal example do? Load the sliced data and check to the corresponding subset of data after full loading of all data?

I am not sure whether this is already the case in RawDataAcquisition, but is the partial loading also possible for other dimensions (diffusion directions, volume repetitions, multiple echoes etc.)?. As you can see in the Data Collection related to our paper, I even had to split a 100 volume ISMRMRD file into 3 smaller ones with 10 GB each (32 coils, 36 slices, long readouts a 60 ms, but nothing unheard of for high-field fMRI). It easily gets unwieldy to keep a reasonable scan in memory, I think, and typically, for the image reconstruction, only 2 to 4 dimensions are used together, whereas the rest is split into different reconstruction tasks.

All the best, Lars

tknopp commented 2 years ago

Ok, so its not just about the slice dimension but it would be helpful to also select partial data for other dimensions. Then we should use keyword arguments that are passed to the RawDataAcquisition constructor. Could you specify which dimensions from the EncodingCounters we want to be able to filter?

Base.@kwdef mutable struct EncodingCounters
  kspace_encode_step_1::UInt16 = 0
  kspace_encode_step_2::UInt16 = 0
  average::UInt16 = 0
  slice::UInt16 = 0
  contrast::UInt16 = 0
  phase::UInt16 = 0
  repetition::UInt16 = 0
  set::UInt16 = 0
  segment::UInt16 = 0
  user::NTuple{8,UInt16} = ntuple(i->UInt16(0),8)
end
tknopp commented 2 years ago

The ugly part is the need to change rawData.params as you already do. I fear this is necessary.

Basically one can use your code as the blueprint for moving that stuff into RawDataAcquisition.

The most crucial part, however, will be to change this line: https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/master/MRIFiles/src/ISMRMRD/ISMRMRD.jl#L21 We need to only read the headers of the profiles. But I am actually not sure if HDF5.jl allows this right now.

Edit: Partial reading seems not be supported by the regular "high-level" HDF5.jl routines but according to this C example it should be doable using the low-level HDF5 API: https://www.asc.ohio-state.edu/wilkins.5/computing/HDF/hdf5tutorial/examples/C/h5_compound.c

alexjaffray commented 2 years ago

The most crucial part, however, will be to change this line: https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/master/src/IO/ISMRMRD.jl#L20 We need to only read the headers of the profiles. But I am actually not sure if HDF5.jl allows this right now.

@tknopp can you update the link? it's broken after #72 was merged.

tknopp commented 2 years ago

done