hMRI-group / hMRI-toolbox

A toolbox for quantitative MRI and in vivo histology using MRI (hMRI)
GNU General Public License v2.0
59 stars 44 forks source link

Toolbox expects SE/STE input sorted by echo and flip angle whereas BIDS sorts by flip angle then echo #46

Closed lukeje closed 1 year ago

lukeje commented 2 years ago

Moving discussion from #45 here.

Another important thing, hMRI-toolbox expect the input of B! maps sorted by echo and flip (in the same order as in raw dicom names): echo-1_flip-01, echo-2_flip-01, echo-1_flip-02, echo-2_flip-02, so be careful when you add them.

P.S. @lukeje The echo/flip sorting is extremely annoying, may be we need to fix it somehow

lukeje commented 2 years ago

If I understand the issue that @nbeliy has correctly, then it is that the expected order of input for the hMRI toolbox is based on an assumed ordering derived from Siemens DICOMs that goes:

  1. set-1_echo-1
  2. set-1_echo-2
  3. set-2_echo-1
  4. set-2_echo-2

whereas the field names specified in BIDS would give the order:

  1. echo-1_flip-01
  2. echo-1_flip-02
  3. echo-2_flip-01
  4. echo-2_flip-02

This seems like it should be fixed by simply sorting the input by echo time (which should always be defined, I think).

nbeliy commented 2 years ago

I think it's an issue of mostly documentation. The hMRI toolbox don't tell explicetly in which order the maps are passed.

Sorting files internally in hMRI toolbox is difficult -- in current implementation the metadata is retrieved only from the first file.

Solutions:

Personally I prefer second solution, because it's explicit and transparent. But, it will require the change of subjob parameter change.

From technical point of view, all three solutions are not complicated to implement.

lukeje commented 2 years ago

Could you check whether 8a4a7de fixes the issue for you? It reads in the echo times for each of the volumes and tries to distinguish "bids order" from "standard order" based on that.

nbeliy commented 2 years ago

It should fix, but it may be not robust agains weird bids names. If you recover the list of echo times, you can sort files by them directly, without assuming the order in which they comes.

Something like:

ET = unique(b1map_params.b1acq.EchoTimes, 'sorted');
V_SE = V(find(b1map_params.b1acq.EchoTimes == ET(1)));
V_TSE = V(find(b1map_params.b1acq.EchoTimes == ET(2)));
lukeje commented 2 years ago

Thanks; that's a good suggestion which I'll try and implement.

lukeje commented 2 years ago

OK, this branch should now fix the problem relatively generally: https://github.com/hMRI-group/hMRI-toolbox/tree/SDAM It'll be a little while till the PR (#41) gets integrated into the main branch (probably next month), but you can already pull the changes into your local repository. I'll leave this issue open for now until the branch is pulled in.

nbeliy commented 2 years ago

Will it be also usefull to retrieve also flip angles from all files, same as echo times?

With failsafe if retrieved FA is 0, then trying to use B1mapNominalFAValues.

The FlipAngle is required in bids, while FlipAngleSeries is custom addition to get full list of fa.

lukeje commented 2 years ago

I wouldn't want to rely on FlipAngle being 0, as this may be entirely specific to a given sequence (e.g. FlipAngle could instead be the excitation flip angle, which is completely different from the "spin echo" flip angle). This again goes back to decisions needing to be made on how to represent this information in the sequence and how to extract this information from the image header files (#16).

nbeliy commented 2 years ago

FlipAngle could instead be the excitation flip angle, which is completely different from the "spin echo" flip angle

Ok then, even if it seems little bit strange for me.

For classic metadata: What I would propose (not in this issue but more general), is to have a separate library in which for each known (or described) sequence have associated values for metadata. get_metadata_val check the file sequence name, cache the info and uses it. Adding new sequence should be then as easy as just add new entry to library.

For bids metadata: They are changing standard again, soon (in their dev branch) the values like FlipAngle or EchoTime can be an array of numbers corresponding to values of all files in that they call file collection. It should simplify a little the retrieval of metadata for b1 maps, but may brake other ones (it will retrieve an array instead of value).

If needed I can try to provide prototypes for these issues.

SiyaMRPhy commented 2 years ago

This issue was giving me trouble. I named echo after flip
sub-BAL001_ses-01_acq-rescanHighFA_flip-01_echo-01_TB1EPI.json sub-BAL001_ses-01_acq-rescanHighFA_flip-01_echo-02_TB1EPI.json sub-BAL001_ses-01_acq-rescanHighFA_flip-02_echo-01_TB1EPI.json sub-BAL001_ses-01_acq-rescanHighFA_flip-02_echo-02_TB1EPI.json

for a flip there are two echos However, it's not a recommended naming. BIDS decided to list entities based on alphabetical order (they arbitrarily decided on this order and its difficult to change it)

One solution is to use https://github.com/bids-standard/bids-matlab

use bids query to grep the correct file and feed to the input of the code in the correct order.

lukeje commented 2 years ago

@SherS2, Does the solution I proposed (which involves sorting by echo time; see PR #41) not work for you?

SiyaMRPhy commented 2 years ago

@lukeje I think that your solution is better than mine.

lukeje commented 2 years ago

@SherS2 Thanks; I hope we will have time to get round to that PR at the PR meeting in August.

SiyaMRPhy commented 2 years ago

@lukeje sure. We can discuss this during PR meeting

I was thinking about using the 'bids query' to feed the input files to the matlabbatch in the correct order.

You are implementing the solution inside the toolbox. That is the ideal solution. I don't need to worry about the order while putting inputs into the toolbox.

nbeliy commented 2 years ago

I was thinking about using the 'bids query' to feed the input files to the matlabbatch in the correct order.

I think it's the best procedure, the bids-matlab queries, but don't forget to crosscheck the order in which bids-matlab returns the paths, also control how much of files you receive. Also, be aware that bids-matlab can be slow to generate database for large datasets.

lukeje commented 1 year ago

@nbeliy @SherS2 Unfortunately I have come across datasets (namely, all the datasets measured at the MPI-CBS in Leipzig) where the lower TE is the stimulated echo, which causes the TE-checking to break. I have thus decided to make it optional for now. I have added an example config file to the branch which shows how to set the flag if you want the TE-checking turned on.

Since we are taking this path for the TEs, would it be worthwhile to do something similar for checking the flip angles in the metadata?

As background, this problem occurs because, in common implementations of the 3D-SE/STE B1 mapping sequence, the TEs that are written into the DICOM headers are those which are chosen in the protocol editor on the scanner. However only the first TE is used to set the timing for the sequence, and so the second TE can be chosen arbitrarily. If the second TE is not manually set to be longer than the first, then sorting by TE will give the "wrong" order. Note that this also seems to be the case for the reference dataset, where the second TE of 130 ms is not in any real sense the TE of the stimulated echo, and so it is just a coincidence that it is longer than the first TE.

A suggestion I would have going forward is that the TE written to the header of the STE DICOMs should be the time since excitation (i.e. the first TE + the mixing time). This would also have the benefit of encoding the mixing time in the header. If I get time I will make this suggestion on the qMRI-BIDS repository...

SiyaMRPhy commented 1 year ago

@luke I use my own script to generate BIDS (in fact BIDS-like, the filename follows BIDS, json is same as that of the original file )

I use acqpar.AcquisitionNumber and acqpar.EchoNumbers to name the files in BIDS format.
example base_name = ['sub-' sub_name ... '_ses-' session ... '_acq-' acq_type ... '_flip-' sprintf( '%02d',l_list.acqpar.AcquisitionNumber)... '_echo-' sprintf( '%02d',l_list.acqpar.EchoNumbers)... '_TB1EPI'];

if acqpar.AcquisitionNumber and acqpar.EchoNumbers are available in the JSON of the BIDS JSON file, then you could sort it easily (also verify the order before processing). @nbeliy might be able to comment on how easy it is to add this in the BIDSme for BIDS conversion

nbeliy commented 1 year ago

@lukeje The ideal solution is to split input files in 2 groups: SE and STE, thus removing the confusion. Unfortunately it will require change in function signature and compatibility issues.

An alternative is to assume (again dangerous) that first echo number encountered is SE, and second is STE (stable instead of sorted):

ET = unique(b1map_params.b1acq.EchoTimes, 'stable');
V_SE = V(find(b1map_params.b1acq.EchoTimes == ET(1)));
V_TSE = V(find(b1map_params.b1acq.EchoTimes == ET(2)));

In both naming schemas, the echo index goes in acquisition order.

In any case, that should be given in documentation (mbatch).

lukeje commented 1 year ago

@SherS2 The aim of the changes is to make it so that you can just give all the qMRI-BIDS compatible SE/STE data as input and the toolbox will sort it and extract the relevant metadata. This should be the case in the latest version of #41 if you modify the b1 processing metadata input in line with this file.

@nbeliy I guess you mean changing the batch input system so that the SE and STE files are separate fields? I don't want to do that as it will (a) require the user to be able to tell which is which and (b) break all current scripts. I have thus gone for the approach mentioned in my response to Siya.

@SherS2 @nbeliy Is there already a public qMRI-BIDS-compatible 3D-EPI SE/STE B1 mapping dataset I can test it out on?

lukeje commented 1 year ago

@SherS2 @nbeliy I found the BIDS-ified reference MPM dataset, so I'll test it out on that.