Support 2D COSY TWIX reader

[darrencl]

Hi, I would like to move our work which parsing 2D COSY Siemens into suspect. What we need is another metadata to record dwell time on T1 axis. Also 2D needs to be able to infer the frequency along both axes. That said, we could make a class for 2D COSY to accommodate those (and other stuff we might find along the way).

After adding that, my plan is to retain the sequence name information in TwixBuilder.

suspect/suspect/io/twix.py

Lines 64 to 72 in 4d2c347

	mrs_data = MRSData(data,
	self.header_params["dt"],
	self.header_params["f0"],
	te=self.header_params["te"],
	tr=self.header_params["tr"],
	metadata=metadata,
	transform=self.header_params["transform"])
	return mrs_data

Then in TwixBuilder.build_mrsdata() above, I am going to put something like:

if self.sequence_name in [".....SVS_SE"]:
    return MRSData(....)
elif self.sequence_name in ["....COSY", ...]:
    return COSYMRSData(...)

This approach allows us to know exactly the sequence, hence, we can be confident if it is a 1-dimensional or not. That said, upon the returning MRSData above, we can correct/remove the unused extra dimension as written in Gannet (also found this issue when loaded on suspect).

See this comment (the dimension I am talking about is referred to the fourth dimension in that comment): https://github.com/richardedden/Gannet3.0/blob/534ce01c824359f84ca48fd17c1562fb2a8f7fd1/SiemensTwixRead.m#L305

With the presence of this COSY class, we might also include DICOM COSY reader, although it might be very dumb compared to this TWIX COSY reader. Since the DICOM saves COSY in individual file per increment, we might need different function that accepts directory (e.g. suspect.io.load_siemens_dicom_2d(dicom_volume)). But of course, let's leave this for later. :-)

Any thoughts?

[bennyrowland]

I am very happy to get COSY better supported in Suspect, and this has come at an interesting time because there are other changes I would like to make relating to axes in general. However, it is important to figure out the best way to do this, which I try to make the most generic way is possible. Does the Twix file for COSY have any information in it to tell us which axis is the indirect frequency axis, or what the time increment is for the evolution time? It is always hard to label the data axes for spectroscopy data because they are often repurposed from the imaging labels for alternative MRS only purposes like editing.

I believe the issue described in Gannet is for the newer Siemens sequences which include prospective frequency drift correction where they periodically collect signal following the water suppression pulse to assess changes in B0 and adjust the TX/RX frequencies to follow it. This data then goes into another axis of the dataset rather than a separate file or measurement, I suspect because it is collected interleaved with the main data.

What I would like to do for Suspect moving forward is to take inspiration from the xarray library to add an Axis object per dimension of the MRSData object. This would be useful in keeping track of e.g. which axis is dynamics and which coils by simply labelling the axes, so you could get a dimension by label and not by index, but it could also be used to provide values for each index along the axis. This would be the time axis for the FID axis, simple channel number for the coil axis etc. It could also be useful for handling FTs, both spectral and spatial, as you could tell whether your X-axis was kx or x to check whether your CSI data was in k-space or image space, for example. So I would imagine that rather than have a COSY class, we would just need to make sure that the appropriate axis got labelled with the evolution time for each axis step. COSY processing functions would just have to check that the relevant axes were available on the supplied MRSData object, probably using a decorator to automate this (I did this for transforms for certain MRSBase functions).

Obviously this is a reasonably big addition to Suspect, although it should be pretty backwards compatible, and there are still some questions about how to implement it, and not least to figure out how to create the relevant axis objects from raw data files like Twix. I would certainly welcome your thoughts on this.

A related project that I have been considering is moving from the current model, where I subclass numpy's ndarray to create MRSBase, to an encapsulation model where MRSBase contains an ndarray as a data object. The disadvantage with this is that you have to maintain a much bigger api, because you have to forward all relevant function calls (like e.g. adding two spectra together) onto the internal object, and so on, but the advantage is that you gain more control over that API, which reduces the chance that your metadata and data get out of sync. For example, if you want to make it so that getting a single voxel out of a CSI dataset will update the transform to reflect the location of just that voxel. I am still trying to decide if this is a good idea or not.

[darrencl]

@bennyrowland As far as I know, TWIX just dumps all the acquisition, but record these information to organize in the loop counters. Hence, we can loop through the acquisition and look into the loop counters to figure out the t1 and average index. @c42f knows more on this since he implements most of the IO stuff in our Julia code. Feel free to jump in if what I said is not right! The dwell time in T1 axis depends on the sequence, for instance SR COSY dumps this in sWipMemBlock.adFree[1], but sometimes this doesn't exists (we just hard code it to 0.8ms when that happens since it is the only value we have seen). As for the number of increments, it can be inferred through repetition in loop counter. Therefore this axis will be incrementing every dwell time T1 from 0.

I agree! xarray seems neat for this purpose. We actually implemented similar thing in Julia using AxisArrays.jl (arrays with axis). I saw xarray library also allows to use the usual slicing (http://xarray.pydata.org/en/stable/indexing.html#indexing-and-selecting-data) by index, so I think there should be no issue in backward compatibility?

I think moving the ndarray as data is a good idea if you're planning to support that example of yours. Although, you're right on making API very big, but I think those functionality just needs to be implemented once and very rarely to be changed? If we're going to change the ndarray to xarray, this is a good time to think about this also.

And yes this will be a quiet big API update. Do you have any thoughts on how we should break this down? I am happy to be involved in this update. :-)

[bennyrowland]

I think this update will need some more thinking about. There is currently an ISMRM working group discussing how to implement a standard data format for MRS data to replace twix/P-files/sdat etc. and I think the decisions reached there will probably mirror what we want to see from an updated MRSData object. @darrencl, if you would like to join this group I am sure you would have some useful contributions, let me know and I will send you an invite.

For example, there is some discussion about sparsely sampled axes (e.g. compressed sensing) where an axis would have to store the values of the relevant parameter for each index, rather than assuming the axis is defined by a first value, a delta and a number of elements. In fact, we can also imagine situations where a single dimension of data defines variations in multiple parameters, for example spiral EPSI where we have a series of k_x/k_y(/k_z?) tuples rather than separate axes for the different spatial directions, or possibly in a multi-spectral-dimensional situation where again we have two (or more) coordinates to be defined for each FID acquired. If the new format doesn't support these situations then we exclude anybody that is using those techniques which kind of missed the point of a standard. On the other hand, in the limit this becomes a 3D array with dimensions of ADC points, coil channels, and a tuple of all loop counter values which puts the burden mainly back on the user again which doesn't help them very much.

I think that using xarray is probably the right choice to add the necessary axis labelling to the ndarray, the question is then whether to do an encapsulation or an extension (http://xarray.pydata.org/en/stable/internals.html#extending-xarray). The encapsulation would be more work but probably provide a nicer interface, and might also allow some extra control if that were necessary.

Probably the next step here is to think about some hypothetical code that we might imagine writing to use these new features, that will guide what the features need to offer. Since you are likely to be the first person to do something different to the existing features, would you like to suggest how you would like to use this COSY object in downstream code?

[c42f]

For example, there is some discussion about sparsely sampled axes (e.g. compressed sensing) where an axis would have to store the values of the relevant parameter for each index, rather than assuming the axis is defined by a first value, a delta and a number of elements. In fact, we can also imagine situations where a single dimension of data defines variations in multiple parameters, for example spiral EPSI where we have a series [...]

My opinion is that a lower-level API should be offered which gives direct access to the sequence of FIDs along with any metadata per FID. This is the only thing which is truly flexible from an MR physics and sequence design perspective.

This is the reason I implemented the MRExperiment type (as a list of Acquisitions) in https://github.com/TRIImaging/MagneticResonanceSignals.jl/blob/1d81be222e247e2b3af030d07924965ec52cb8bb/src/io_twix.jl#L162

If you look at the Julia code, you'll see that the higher level interpretation of the signal as a cartesian sampling of the delay time T1 and phase cycling dimensions in COSY are just a lightweight view on top of the list of acquisitions.

The problem, as you've pointed out, is that there's an endless set of design possibilities for sequences. And I think trying to fully categorize these is ultimately a hopeless task. Of course, categorizing a subset and providing a more convenient API is fine and good! But the lower level API should exist for true flexibility.

(In my opinion we've only seen the start of weird hard-to-categorize sequences with things like MRF. In the future (or maybe already; I haven't been following the research) I think we'll see acquisition and reconstruction become much more tightly coupled with on-the-fly sequence adjustment to maximize the task-specific information gained per unit time. Presumably using some composition of physics-based and data driven models. This is the research I'd want to pursue if I ever got back into the field :-) )

[darrencl]

@bennyrowland I'm so sorry I missed this thread! Have you made any decision on how should we approach this?

I am happy to join that group, although I am not sure if I can be of any help, since i am pretty much only working with single voxel spectro data. Feel free to send me an invite! :-)

But the lower level API should exist for true flexibility.

+1 on this!

[bennyrowland]

@darrencl sorry for taking so long to respond on this thread - I have been focussed on other things a long way from spectroscopy for a while. I have however been playing a lot with xarray as a potential backend for spectroscopy data - it gives a lot of power, particularly with labelling dimensions, but also has some challenges around preserving metadata (TE/TR, transform etc.) I think this is the right way to go in the longer term, but it will break a bunch of stuff for a while. I am hoping to get some time to work on this in a while, will probably be in a branch where I can play with stuff a bit more freely.