MRIBase.rawdata does not work correctly with shuffled profiles

[atsanda]

In the current form MRIBase.rawdata operates under assumption that profiles come in a certain order which then allows to place readout data in certain places of the multidimensional kdata array. I find this assumption relatively strong and suggest lifting it. I am going to write a test suite to highlight the issue and to later fix the bug.

[atsanda]

Here is an example:
https://github.com/MagneticResonanceImaging/MRIReco.jl/blob/bug-issue-215-MRIBase.rawdata-shuffled-profiles/MRIBase/test/testRawAcqData.jl

[aTrotier]

Do you expects MRIBase.rawdata to return something ordered like a k-space ?

This is done at the next steps after the conversion from raw -> acq and then you can use kDataCart(acq::AcquisitionData)

[tknopp]

yes, rawdata should basically output the data in the way that it has been acquired by the MRI scanner. The resorting is done afterwards when converting to an AcquisitionData object.

[atsanda]

@aTrotier, yes, kspace should be ordered after that, I mean the encoded part of it. I've taken a look into kDataCart. I don't see how it could help solving the issue. It orders according to repetitions, slices, echos and channels but the data within profiles is just copied. Do I understand it right?

[aTrotier]

@atsanda MRIBase.rawdata is not suppose to do that.
The RawAcquisitionData structure is mostly a copy of the ISMRMRD standard which stores the data in the order they are acquired. The rawdata function just retrieve the rawdata for specific echo, slice and rep which right now are the 3 dimensions others than x/y/z which are used later by MRIReco

[atsanda]

@aTrotier thank you for pointing me to the kDataCart. I've tested it and it fixes the order nicely. I've put our discussion into test cases for RawAcqData and AcqData. I am going to create a pull request and close the issue as "not a bug".

Yet, the reason I started to blame the rawdata function lies in its complexity. I don't think the function follows the guidelines mentioned by you and Tobias. Namely, it does filtering unique profiles, distinguishes between "custom" and other trajectories, allocating dense array differently for the latter one. All of that made me think that this function can be "smarter" than it is. If I had to implement it on my own based on your description I'd end up with something like this:

function kspace_idx(f::RawAcquisitionData, slice, contrast, repetition)
  sl = slices(f)
  rep = repetitions(f)
  contr = contrasts(f)
  idx_sl = findall(x->x==slice,sl)
  idx_contr = findall(x->x==contrast,contr)
  idx_rep = findall(x->x==repetition,rep)
  idx = intersect(idx_sl,idx_contr,idx_rep)
  return idx
end

function rawdata(f::RawAcquisitionData; slice, contrast, repetition)
  idx = kspace_idx(f, slice, contrast, repetition)

  # assume the same number of samples for all profiles
  numSampPerProfile, numChan = size(f.profiles[idx[1]].data)
  numSampPerProfile -= (f.profiles[idx[1]].head.discard_pre+f.profiles[idx[1]].head.discard_post)

  kdata = zeros(typeof(f.profiles[1].data[1, 1]), numSampPerProfile, length(idx), numChan)
  for (i, profile) in enumerate(f.profiles[idx]) 
    i1 = proile.head.discard_pre + 1
    i2 = i1 + numSampPerProfile - 1
    kdata[:, i, :] .= profile.data[i1:i2, :]
    if flag_is_set(profile, "ACQ_IS_REVERSE")
      kdata[:,i,:] .= reverse(kdata[:,i,:], dims=1)
      flag_remove!(profile,"ACQ_IS_REVERSE")
    end
  end

  return reshape(kdata, :, numChan)
end

Which would always allocate the number of profiles found for a certain slice, rep, contrast w/o doing smart calulations or filtering. And since I've spent some time understanding this function, my question to you would be if we should refactor it this way making it nice and clean? Unfortunately, the function is not covered with tests to check why this complexity appeared there in the first place.

And I also have another point. I started looking into that for a reason: I have 3d MRI data where profiles have weird order + along y and z kspace is not rectangular but has ellipsoidal shape. And I try to apply direct reco to that using reconstruct function which fails. kDataCart indeed helps but it is not called. Instead kData is called here:

            kdata = arrayType(kData(acqData,k,j,i,rep=l)) .* arrayType((weights[k].^2))

Then it is just the output of the rawdata which is taken for the reconstruction. So, in this case, should I manually extract kspace with kDataCart and pass it to inverse fourier or is there a standard way of doing it?

[aTrotier]

kDataCart is a function I wrote first to help me for conversion back and forth in order to call BartIO.jl
I think it is only used in espirit for simplicity.

kData is really just extracting the kdata for a certain k,j,i,rep

"""
    kData(acqData::AcquisitionData, echo::Int64=1, coil::Int64=1, slice::Int64=1;rep::Int64=1)

returns the k-space contained in `acqData` for given `echo`, `coil`, `slice` and `rep`(etition).
"""
function kData(acqData::AcquisitionData, echo::Int64=1, coil::Int64=1, slice::Int64=1;rep::Int64=1)
  return acqData.kdata[echo,slice,rep][:,coil]
end

but the encoding operator F is used here to perform the FFT and subsample / put the corresponding k-space lines at the right place. The ordered k-space is not really created anywhere.

You're right the rawdata function is especially complicated but I am not sure we can get rid of the filter. We should add a test to see what happens if a line is present 2 times in the profiles.

I guess in that case the data should be averages but in the current implementation we only use the first one

[aTrotier]

Then it is just the output of the rawdata which is taken for the reconstruction. So, in this case, should I manually extract kspace with kDataCart and pass it to inverse fourier or is there a standard way of doing it?

      if !MRIBase.isUndersampledCartTrajectory(shape,tr)
        ftOp = FFTOp(Complex{T}; shape, unitary=false, S = S, fftParams(S)...)
      else
        idx = MRIBase.cartesianSubsamplingIdx(shape,tr)
        ftOp = SamplingOp(Complex{T}; pattern=idx, shape, S = S) ∘ FFTOp(Complex{T}; shape, unitary=false, S = S, fftParams(S)...)
      end

@atsanda Humm, actually if the k-space is not detected as UnderSampledCartTrajectory I suppose it is a simple FFT operator

[atsanda]

@aTrotier

You're right the rawdata function is especially complicated but I am not sure we can get rid of the filter. We should add a test to see what happens if a line is present 2 times in the profiles.

I would expect it to fail due to dimensions mismatch with the Fourier operator. Data would have more points then than the operator would expect. I think the operator is allocated according to the trajectory which would contain the shape of the kspace. It is rather the question of responsibility for this function. We can keep it simple and focused moving filters to the point where AcqData is created and where we expect the data to be transformed. Then it should work. I can prepare a pull request for that with tests if you agree that it is a good idea. I just need a point somewhere in the code where constructing a proper dense array for kspace with filters and actual indexes would be expected. If that's not the rawdata, would be great to find one. Otherwise the framework does not work for me.

Humm, actually if the k-space is not detected as UnderSampledCartTrajectory I suppose it is a simple FFT operator

And that is my case, I just don't have some points around the y-z corners in my ksapce and I am fine with zeros there.

[aTrotier]

function isUndersampledCartTrajectory(shape::NTuple{3,Int}, tr::Trajectory)
  return shape[1]!=numSamplingPerProfile(tr) || shape[2]!=numProfiles(tr) || shape[3]!=numSlices(tr)
end

in your case the function returns true. can you remove the if part and always apply the samplingOp:

      if !MRIBase.isUndersampledCartTrajectory(shape,tr)
        ftOp = FFTOp(Complex{T}; shape, unitary=false, S = S, fftParams(S)...)
      else
        idx = MRIBase.cartesianSubsamplingIdx(shape,tr)
        ftOp = SamplingOp(Complex{T}; pattern=idx, shape, S = S) ∘ FFTOp(Complex{T}; shape, unitary=false, S = S, fftParams(S)...)
      end

by

        idx = MRIBase.cartesianSubsamplingIdx(shape,tr)
        ftOp = SamplingOp(Complex{T}; pattern=idx, shape, S = S) ∘ FFTOp(Complex{T}; shape, unitary=false, S = S, fftParams(S)...)

@tknopp Should we force the FFTOp to always be multiply by the SamplingOp ? Or do you want to keep the if/else and change how is define an undersampledCartTrajectory by something like :

idx = MRIBase.cartesianSubsamplingIdx(shape,tr) 
isUnderSample = length(idx) < prod(shape) ? true : false

[tknopp]

I would not always perform the subsampling. It is an unnecessary operation for fully sampled data.

[aTrotier]

the conversion to AcquisitionData already filter the supplementary profiles (for cartesian acquisition at least) :

using MRIReco
include("MRIBase/test/Utils.jl")

n_profiles = 10
n_samples_per_profile = 12
data = rand(ComplexF32, n_samples_per_profile, n_profiles)
f = Utils.get_default_cartesian_raw_acq_data_2d(n_profiles, n_samples_per_profile, data)

acq_1 = AcquisitionData(f)
length(acq_1.kdata[1,1,1]) == n_profiles * n_samples_per_profile

f2 = deepcopy(f)
push!(f2.profiles,f2.profiles[1])
acq_2 = AcquisitionData(f2)
length(acq_2.kdata[1,1,1]) == n_profiles * n_samples_per_profile
kDataCart(acq_2) == kDataCart(acq_1) # equal because only the first profile is used

f3 = deepcopy(f)
prof = deepcopy(f3.profiles[1])
prof.data .= 1
push!(f3.profiles,prof)
acq_3 = AcquisitionData(f3)
length(acq_3.kdata[1,1,1]) == n_profiles * n_samples_per_profile
kDataCart(acq_3) == kDataCart(acq_1) # equal because only the first profile is used 
kDataCart(acq_3) == kDataCart(acq_2) # equal because only the first profile is used 


f4 = deepcopy(f)
push!(f4.profiles,f4.profiles[1])
f4.profiles[1].data .= 1
acq_4 = AcquisitionData(f4)
length(acq_4.kdata[1,1,1]) == n_profiles * n_samples_per_profile
kDataCart(acq_4) == kDataCart(acq_1) # false because the first profile is used and is equal to 1

[atsanda]

Hm, sounds like a good reason to remove unnecessary filtration then. I will look into it and create a separate issue for refactoring.

Regarding the subsampling operator, it can help me with missing profiles but does not solve the issue with the ordering which I wanted to solve with this issue. As kData just takes the output of rawdata it is assumed that the order of profiles matches the memory layout of multidimensional kspace. It is not the case for me. kDataCart fixes that with proper re-indexing. I think we should integrate that into the creation of AcqData struct. it should not hurt the package as it makes it more general. That would be the second refactoring issue. @aTrotier , what would you say to that?

[aTrotier]

I still don't understand what is the issue with the indexing.

Yes, the kData are not reindexed (they are not supposed to). Even if we reordered them during the conversion to AcquisitionData, if the data are undersampled (with a CS or IPAT acceleration for example), the "order" kdata won't be useful. It will works only for a fully sampled kspace without acceleration, zero-filling, partial-fourier, cropped circle encoding.

[atsanda]

Well, first of all, this concrete issue, as discussed above, should be closed. MRIBase.rawdata is not supposed to do anything with profiles (although it does in fact). I wrote tests to document that. @aTrotier, could you, please, review the PR and merge? We could then close here with the rawdata.

However, the issue I have is still there. I have a 3D Cartesian measurement where kpace was not measured around the corners of ky-kz, a sort of ellipsoidal mask was applied along ky-kz. I try to reconstruct it with the framework and it fails. It happens because direct reconstruction accesses rawdata to apply inverse Fourier. The rawdata recognizes Cartesian trajectory, pre-allocates a full zero-filled kspace for that and fills that in the order profiles appear. In my case, that results in a complete corrupted kspace data. I expect missing parts to be filled with zeros and placed where they belong. I've boiled down this big issue to the problem of ordering. If at some point rawdata or a function afterwards did not just take profiles the way they appear in the rawdata but took into account actual ky and kz values, my issue would be solved. It requires lifting the assumption that profiles come in the right order and there are enough profiles to fill a Cartesian grid.
I can also craft a toy ISMRMRD file to showcase my struggle. It well may be that the issue can be solved differently,