Get correct shape of k-space trajectory tensor

[schuenke]

Currently, the KTrajectory ABC (_KTrajectory.py) has a method _get_shape() to calculate the correct shape of the traj torch.tensor:

class KTrajectory(ABC):
    """Base class for k-space trajectories."""

    @staticmethod
    def _get_shape(header: KHeader) -> tuple[int, ...]:
        """Get the shape of the trajectory for the given header."""
        limits = header.encoding_limits
        other_dim = np.prod(
            [
                getattr(limits, field.name).length
                for field in dataclasses.fields(limits)
                if field.name not in ('k0', 'k1', 'k2', 'segment')
            ]
        )
        shape = (
            other_dim,
            3,
            limits.k2.length,
            limits.k1.length,
            limits.k0.length,
        )
        return shape

A discussion came up if this still works for all edge cases like acquisitions with freq/phase oversampling, partial fourier etc.

I did some test measurements with a simple FLASH sequence that @hdillinger compiled with the following default settings:

• FOV 192 x 192 x 5 mm³
• matrix size 192 x 192 x 1

I converted the Siemens XA61 raw data files to ISMRMRD files using the latest siemens_to_ismrmrd converter, read in the data using the MRPro KData class and extracted the following infos from the corresponding KHeader objects:

for header in header_list:
    print(f'name: {header.protocol_name} - enc matrix: {header.encoding_matrix} - recon matrix: {header.recon_matrix}')

name: FLASH_basic - enc matrix: SpatialDimension(x=384, y=192, z=1) - recon matrix: SpatialDimension(x=192, y=192, z=1)
name: FLASH_phase_oversampling_20 - enc matrix: SpatialDimension(x=384, y=232, z=1) - recon matrix: SpatialDimension(x=192, y=192, z=1)
name: FLASH_phase_resolution_75 - enc matrix: SpatialDimension(x=384, y=192, z=1) - recon matrix: SpatialDimension(x=192, y=192, z=1)
name: FLASH_phase_partialfourier_6_8 - enc matrix: SpatialDimension(x=384, y=192, z=1) - recon matrix: SpatialDimension(x=192, y=192, z=1)
name: FLASH_avg4 - enc matrix: SpatialDimension(x=384, y=192, z=1) - recon matrix: SpatialDimension(x=192, y=192, z=1)
name: FLASH_slices3 - enc matrix: SpatialDimension(x=384, y=192, z=1) - recon matrix: SpatialDimension(x=192, y=192, z=1)
name: FLASH_measures5 - enc matrix: SpatialDimension(x=384, y=192, z=1) - recon matrix: SpatialDimension(x=192, y=192, z=1)
name: FLASH_GRAPPA2 - enc matrix: SpatialDimension(x=384, y=192, z=1) - recon matrix: SpatialDimension(x=192, y=192, z=1)

for header in header_list:
    print(f'name: {header.protocol_name} - {header.encoding_limits.k1} (length = {header.encoding_limits.k1.length})')

name: FLASH_basic - Limits(min=0, max=191, center=96) (length = 192)
name: FLASH_phase_oversampling_20 - Limits(min=0, max=230, center=115) (length = 231)
name: FLASH_phase_resolution_75 - Limits(min=0, max=143, center=72) (length = 144)
name: FLASH_phase_partialfourier_6_8 - Limits(min=0, max=143, center=96) (length = 144)
name: FLASH_avg4 - Limits(min=0, max=191, center=96) (length = 192)
name: FLASH_slices3 - Limits(min=0, max=191, center=96) (length = 192)
name: FLASH_measures5 - Limits(min=0, max=191, center=96) (length = 192)
name: FLASH_GRAPPA2 - Limits(min=0, max=190, center=96) (length = 191)

for kdata in kdata_list:
    print(f'name: {kdata.header.protocol_name} - "_get_shape" return: {DummyTrajectory()._get_shape(kdata.header)}')  

name: FLASH_basic - "_get_shape" return: (1, 3, 1, 192, 384)
name: FLASH_phase_oversampling_20 - "_get_shape" return: (1, 3, 1, 231, 384)
name: FLASH_phase_resolution_75 - "_get_shape" return: (1, 3, 1, 144, 384)
name: FLASH_phase_partialfourier_6_8 - "_get_shape" return: (1, 3, 1, 144, 384)
name: FLASH_avg4 - "_get_shape" return: (4, 3, 1, 192, 384)
name: FLASH_slices3 - "_get_shape" return: (3, 3, 1, 192, 384)
name: FLASH_measures5 - "_get_shape" return: (5, 3, 1, 192, 384)
name: FLASH_GRAPPA2 - "_get_shape" return: (1, 3, 1, 191, 384)

Some obvervations/quenstions:

1. it looks like a frequency oversampling by a factor of 2 is set by default. I couldn't find the corresponding field in the sequence card(s) at the scanner. This is most probably the reason for the 384 values that appear in enc matrix.
2. I don't really understand why the recon matrix for phase resolution = 75% and partial fourier = 6/8 is the same. A reduced phase resolution should lead to reduced matrix sizes, no?
3. the k0 dimension of ktraj should NOT include the oversampling data, or?
4. the calculation of "other dim" seems to work for all tested cases

If you want to check the header information on your own, here is a code snippet to create the lists I used above (you might need to disable the shape checks in KData.from_file()):

folder = Path(R'Z:\_allgemein\projects\pulseq\measurements\2023-08-15_T2MES_MRPro\FLASH')
mrd_list = list(folder.glob('*.mrd'))

kdata_list = []
for mrd in mrd_list:
    kdata_list.append(KData.from_file(filename=mrd, ktrajectory=DummyTrajectory()))

header_list = []
for kdata in kdata_list:
    header_list.append(kdata.header)

[ckolbPTB]

Here is some more information regarding the relationship between recon matrix size and encoding limits: https://ismrmrd.readthedocs.io/en/latest/mrd_header.html

Some thoughts:

In my opinion the k-space trajectory should be the exact same size as the k-space data. The shape of the trajectory can therefore either be obtained from kdata.data or from the AcqIdx.k0/k1/k2 tensors.

The encoding limits just tell the min, max and centre value but there is no garantue that all the values between min and max are acquired. One example I think is the GRAPPA2 dataset above. According to the limits there are 191 phase encoding points but in reality there should be less, otherwise it would not be undersampled.

In theory there is nothing which prevents the user from obtaining average 0, 2, 4 and 10. If this makes sense or not is a different question. It would lead to a valid ismrmrd file but our calculation of other_dim would fail.

As far as I remember the oversampling along the frequency encoding direction is always 2 unless the FOV is very large. Then I think the oversampling gets reduced to achieve a max FOV value. The oversampling factor should be the ratio of the encoded matrix size to recon matrix size for the frequency encoding dimension.

The k0 dimension of ktraj should include the oversampled data to ensure kdata and ktraj match. Only once we have removed the oversampled data should the trajectory be adapted.

[fzimmermann89]

Here is some more information regarding the relationship between recon matrix size and encoding limits: https://ismrmrd.readthedocs.io/en/latest/mrd_header.html

Also, the ISMRMRD paper ( https://onlinelibrary.wiley.com/doi/10.1002/mrm.26089) has some information, examples and figures. @SherineBrahma volunteered to provide a short overview of the meanings of the different trajectory-related fields in one of the next meetings.

In my opinion the k-space trajectory should be the exact same size as the k-space data.

This is what we currently enforce and should keep enforced in the KData class, yes.

The shape of the trajectory can therefore either be obtained from kdata.data or from the AcqIdx.k0/k1/k2 tensors.

I would prefer to read this information ab initio from the header 'how it should be' and then match it to the kdata.shape. I feel this might be a bit clearer to understand later on and to find issues. What are your opinion @schuenke ?

The encoding limits just tell the min, max and centre value but there is no garantue that all the values between min and max are acquired.

Exactly! I should have removed the _get_shape in the github at all, it just caused confusion an was just meant as a simple way to get the dummy trajectory to work for some initial tests.

But we will need the limits, the acq idx and one of the matrix sizes to figure out the true dimensions and offsets of the Trajectory, yes.

The k0 dimension of ktraj should include the oversampled data to ensure kdata and ktraj match. Only once we have removed the oversampled data should the trajectory be adapted.

Yes. Fully agree.

[schuenke]

Yeah, I agree that we should try to extract the shape info from the header and than enforce that it matches the kdata shape. I think I have a solution that works for all cases. I will post an update soon.

[schuenke]

Okay, I tried to extract the shape from the header infos and came up with the following solutions:

all parameters except k0 should be extracable from acq_info.idx / AcqIdx dataclass using

val = np.unique(header.acq_info.idx.<param>).size

However, for user there are always 2 unique values: 0 and 96 in all examples above.

The k0 dim should be given by the largest readout, no? If so, the dim can be extracted using:

np.amax(header.acq_info.number_of_samples.detach().cpu().numpy(), axis=(0, 1, 2))

or the torch.amax() function.

The results for the example cases (when ignoring "user") are:

name: FLASH_basic - shape(header): (1, 3, 1, 192, 384) - data.shape: torch.Size([1, 34, 1, 192, 384])
name: FLASH_phase_oversampling_20 - shape(header): (1, 3, 1, 231, 384) - data.shape: torch.Size([1, 34, 1, 231, 384])
name: FLASH_phase_resolution_75 - shape(header): (1, 3, 1, 144, 384) - data.shape: torch.Size([1, 34, 1, 144, 384])
name: FLASH_phase_partialfourier_6_8 - shape(header): (1, 3, 1, 144, 384) - data.shape: torch.Size([1, 34, 1, 144, 384])
name: FLASH_avg4 - shape(header): (4, 3, 1, 192, 384) - data.shape: torch.Size([4, 34, 1, 192, 384])
name: FLASH_slices3 - shape(header): (3, 3, 1, 192, 384) - data.shape: torch.Size([3, 34, 1, 192, 384])
name: FLASH_measures5 - shape(header): (5, 3, 1, 192, 384) - data.shape: torch.Size([5, 34, 1, 192, 384])
name: FLASH_GRAPPA2 - shape(header): (1, 3, 1, 108, 384) - data.shape: torch.Size([1, 34, 1, 108, 384])

This is what we expect, no?

[fzimmermann89]

This is different from what I meant...

The information which data is acquired is also supposed to be described in the header without having to resort to counting the acquisitions (i.e. without val = np.unique(header.acq_info.idx.).size) (See the ismrmrd paper)

And I'd still strongly prefer to do this 'properly' and parsing the information in the header instead of fitting our trajectory calculation to the acquisitions/data. (Except maybe for k0/number_of_samples) ;)

Am Mi., 16. Aug. 2023 um 12:58 Uhr schrieb Patrick Schuenke < @.***>:

Okay, I tried to extract the shape from the header infos and came up with the following solutions:

all parameters except k0 should be extracable from acq_info.idx / AcqIdx dataclass using

val = np.unique(header.acq_info.idx.).size

However, for user there are always 2 unique values: 0 and 96 in all examples above.

The k0 dim should be given by the largest readout, no? If so, the dim can be extracted using:

np.amax(header.acq_info.number_of_samples.detach().cpu().numpy(), axis=(0, 1, 2))

or the torch.amax() function.

The results for the example cases (when ignoring "user") are:

name: FLASH_basic - shape(header): (1, 3, 1, 192, 384) - data.shape: torch.Size([1, 34, 1, 192, 384])name: FLASH_phase_oversampling_20 - shape(header): (1, 3, 1, 231, 384) - data.shape: torch.Size([1, 34, 1, 231, 384])name: FLASH_phase_resolution_75 - shape(header): (1, 3, 1, 144, 384) - data.shape: torch.Size([1, 34, 1, 144, 384])name: FLASH_phase_partialfourier_6_8 - shape(header): (1, 3, 1, 144, 384) - data.shape: torch.Size([1, 34, 1, 144, 384])name: FLASH_avg4 - shape(header): (4, 3, 1, 192, 384) - data.shape: torch.Size([4, 34, 1, 192, 384])name: FLASH_slices3 - shape(header): (3, 3, 1, 192, 384) - data.shape: torch.Size([3, 34, 1, 192, 384])name: FLASH_measures5 - shape(header): (5, 3, 1, 192, 384) - data.shape: torch.Size([5, 34, 1, 192, 384])name: FLASH_GRAPPA2 - shape(header): (1, 3, 1, 108, 384) - data.shape: torch.Size([1, 34, 1, 108, 384])

This is what we expect, no?

— Reply to this email directly, view it on GitHub https://github.com/ckolbPTB/mrpro/issues/42#issuecomment-1680391640, or unsubscribe https://github.com/notifications/unsubscribe-auth/ABH3DV3E5DAZUZE5EEMVRK3XVSRVLANCNFSM6AAAAAA3SFVJ4E . You are receiving this because you commented.Message ID: @.***>

[ckolbPTB]

The information which data is acquired is also supposed to be described in the header without having to resort to counting the acquisitions (i.e. without val = np.unique(header.acq_info.idx.).size) (See the ismrmrd paper)

I am not so sure about this and the ISMRMRD paper only discusses a very simple case. E.g. for the GRAPPA2 case we would need to replicate exactly how Siemens calculates the phase encoding points for kernel calibration and the undersampled lines. All of this has then again be adapted by the number of lines in each GRE shot. I think it is easier to utilise AcqInfo.

[fzimmermann89]

Maybe I am misunderstanding something and we can talk about it in the next meeting, but I am still convinced that using the Idx Counters is not the way to go. This will only hide inconsistencies.

If we need more information than provided in the header (i.e. undersampling scheme), this is something we should provide to the Trajektory class and not try to estimate by to number of points we need to calculate to match the number of unique acquisition counters / k space steps.

So using AcqInfo and Header -> Yes. Counting something in AcqInfo.Idx -> No.

Either way, we need to match exactly the number of ACS lines and the phase encoding points for undersamploing! The AcqInfo.Idx counters will not help for that. It will only tell us how many lines there are, but we still have to figure out where they are supposed to be.

So in my head the issue is exactly the opposite: If we use the counters to figure out where the points are, we are - implicitly - fixing ourselves exactly to how Siemens does it. Instead, we should either ask the user for information OR try to use the information provided in the header.

And in the end, the number of calculad points should then match (what we sanity check for) the number of k space data points and thereby also the number of unique acqinfo.idx values.

[ckolbPTB]

So using AcqInfo and Header -> Yes. Counting something in AcqInfo.Idx -> No.

Agreed!

The AcqInfo.Idx counters will not help for that. It will only tell us how many lines there are, but we still have to figure out where they are supposed to be.

AcqInfo.Idx will tell you the exat position of each readout in k-space and AcqIdx.Flags will tell you if it is a ACS line or not. There should not be any need for additional info from the user.

Maybe a more basic question - do we really need to explicitly calculate the shape of the trajectory? Should it not simply be the result of operations on AcqInfo.Idx?

[fzimmermann89]

Yes!!

Sorry if that was not clear. I am arguing for not trying to figure out the 'shape' but to figure out each value of ktraj. And then in the Kdata class we check for the shape to match the data (and thereby the idx)

So it sounds like in principle we are on the same page.

Once again, I should have put get_shape into the dummy trajectory. It has only any use for returning an empty trajectory. I will move it and create a PR to avoid further confusion.

[fzimmermann89]

(my phone's battery died, part 2.. )

AcqInfo.Idx will tell you the exact position of each readout in k-space and AcqIdx.Flags will tell you if it is an ACS line or not. There should not be any need for additional info from the user.

Maybe I am misunderstanding either you or the standard, but my understanding was, that AcqInfo.Idx are indices into the trajectory (not "the exact position"). So by itself, they don't tell you anything but where to look in the Traj data.

So consider 2D Cartesian, phase undersampling 4. The Idx Counter will still count (0,1,2,3,4......) -> not skipping any lines And we still need to know if the first acquired line is at position 0, 1, 2, or 3 in the matrix. And if the pattern is reset after the ACS lines to be symmetric about the ACS, or to be symmetric about the DC, or just continues. I am not sure where we can get this information from the ismrmrd file or if this is something we have to cover in the trajectory classes -- either by having different classes or by having some __init__ parameters.

PS:

AcqIdx.Flags will tell you if it is a ACS line or not

I didn't know that! Good to know!

[ckolbPTB]

Maybe I am misunderstanding either you or the standard, but my understanding was, that AcqInfo.Idx are indices into the trajectory (not "the exact position"). So by itself, they don't tell you anything but where to look in the Traj data.

So consider 2D Cartesian, phase undersampling 4. The Idx Counter will still count (0,1,2,3,4......) -> not skipping any lines

Yes and no. You are right in the sense that Idx Counter is - as the name suggests - an index and hence of integer value. Therefore, e.g. for a radial acquisition it does not tell you the k-space position but only the index of the radial line and the user needs to provide the angle increment to actually calculate the k-space position. Nevertheless, for phase undersampling of 4 the index (at least for Siemens) would actually be 0,4,8,... and hence sufficient for reconstruction.

[ckolbPTB]

Results of final group discussion: Trajectory should fit to data size and no _get_shape() is required.