MR image data representation as NiftiImages

[johannesmayer]

There is a fundamental issue in our approach of trying to convert MR data into the Nifti images objects.

MR data have many different dimensions, obviously 3 spatial ones, but the additional ones (at least defined in ISMRMRD) are:

• contrast
• phase
• repetition
• set

Currently, we try to take a set of MR reconstructions and convert it into a 3D image. This can't work.

When we are trying to convert these into a 3D Nifti image this requires that
a) the images are consistent wrt their content, (e.g. no images of two slices from contrasts should appear in the same image)
2) the multi-slice 2D acquisitions are placed correctly (i.e. slice thickness and slice distance are taken into account separately and correctly)

The good news is: our reconstructions stored in our image containers are ISMRMRD::Image objects so all our reconstructions already contain the correct information in their header. Our attempt of converting it into a single 3D space is, however, flawed.

[evgueni-ovtchinnikov]

We discussed the issue of the data dimensions many times (see e.g. #315) but could never agree on it.
I for one cannot be of much help here, having no background in MR (no idea what contrast/phase/set is), and I could never stomach ideas to work with 11 or so dimensions. I can only use common sense, which tells me that, after all, what we expect from a scanner are images of a 3D object, do we not?

[DANAJK]

The patient is a 3D object, but they can change in time (for example beating heart), or you can inject a contrast agent that changes the image intensities, or you might look at more than one tissue property such as T1, T2, water, fat etc. or you might acquire data multiple times for future averaging.

We should use the power of the computer language to handle the mapping from dimensions to storage in memory. If we try to put everything into a 3D array, we will end up having to do ourselves a lot of conversions of indices from our physical dimensions to this 3D array.

I see broadly two options:
We pre-specify that at certain points in the recon, arrays will have a fixed number of dimensions. The physical meaning of the dimensions will depend on the application and has to be documented. That fixed number will usually be more than 3.
We provide a function that returns the physical meaning of each of the dimensions in the array, and that array is ND with N able to change.

[KrisThielemans]

As discovered in #1186 we currently have a crash when constructing NiftiImageData from MRImageData that is not simply 3D.

At present, conversion to NiftiImageData happens via NiftiImageData::operator=(const ImageData&) . This first creates an empty image

SIRF/src/Registration/cReg/NiftiImageData.cpp

Line 101 in d034e4c

this->_nifti_image = NiftiImageData<float>::create_from_geom_info(*to_copy.get_geom_info_sptr());

and then copies data across. However, the creation is based only on the geom_info which contains no information on repetitions, contrast etc, therefore creating an image with only 1 repetition/contrast etc. The copy will then segfault as it writes in non-allocated memory. #1226 is working around this by only copying data from the first volume, but this is of course not-ideal.

As NiftiImageData can support multi-dim data (up to 7 dimensions I believe), we should try and preserve it. However, Nifti does not have any definitions corresponding to MR data.

Some pointers:

• From the Nifti spec chapter 2:
  • Nifti dim[4] corresponds to time
  • Nifti dim[5] corresponds to "vectors". Meaning is indicated via intent_code etc, see Nifti1 doc on possible values. It is used by many, including NiftyReg and hence SIRF to store deformation/displacements fields. I believe some of this usage is NiftyReg specific (as it goes beyond the standard)
  • I see no mention of dim[6] and dim[7]. Presumably different programs have different conventions.
• We have VoxelisedGeometricalInfo3D which contains info on number of vixels, their size and orientation. Extending this class to contain more info would be possibly but it would need to change name. Probably better to derive class from it with the extra info, or have a MetaDataInfo class that contains a VoxelisedGeometricalInfo3D
• Current implementation of NiftiImageData::create_from_geom_info uses the is_tensor arg to force dim[5] to the size of the deformation vector (2 or 3). The intent is then set in the appropriate constructors, e.g.

  SIRF/src/Registration/cReg/include/sirf/Reg/NiftiImageData3DDisplacement.h

  Lines 72 to 73 in d034e4c

  	NiftiImageData3DDisplacement(const inputType * const data, const VoxelisedGeometricalInfo3D &geom)
  	: NiftiImageData3DTensor<dataType>(data, geom) { this->_nifti_image->intent_code = NIFTI_INTENT_VECTOR; this->_nifti_image->intent_p1=1; }

From this, it seems that it might not be too hard to "collapse" all MR dimensions into the Nifti dynamic field, or use dim 4, 6 and 7. This could be done by something like this

NiftiImageData<dataType>& NiftiImageData<dataType>::operator=(const ImageData& to_copy)
{
 ...
       auto nii_ptr = dynamic_cast<const NiftiImageData<dataType> * const >(&to_copy);
        if (nii_ptr) {
            // Check the image is copyable
            if (!nii_ptr->is_initialised())
                throw std::runtime_error("Trying to copy an uninitialised image.");

            copy_nifti_image(_nifti_image,nii_ptr->_nifti_image);
            this->_data = static_cast<float*>(_nifti_image->data);
            this->_original_datatype = nii_ptr->_original_datatype;
            set_up_geom_info();
        }
        else {
            int num_time_frames = 1;
            int num_repetitions = 1;
            ...
            if (auto mr_ptr = dynamic_cast<const MRImageData<dataType> * const >(&to_copy))
            {
               num_time_frames = ...; etc
            }
           // temp class to store repetition
            VolumeInfo vol_info(num_time_frames, num_repetitions, ...);
            this->_nifti_image = NiftiImageData<float>::create_from_geom_info(*to_copy.get_geom_info_sptr(), vol_info);
           }
}

[KrisThielemans]

@ashgillman maybe you have some more experience with this.

[ashgillman]

Just came across this.

I did read this page recently that I don't think you linked, the documentation on dim: https://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/dim.html

My 2c, I think Kris' suggestion is quite sane: to "ravel"/"unravel" (in numpy language) these different dimensions into one. But not dim 4.
Dim 5 would make sense here, with intent_code=1011, and encoding the ravel key in intent_name. You only get 16 chars, but might be enough with a code for slowest-to-fastest like "SIRF 4c3p8r4s" with "SIRF" as a magic code and 4c3p8r4s meaning 4channels, 3 phase, 8 repeats, 4 sets from slowest-to-fastest varying.

 #define NIFTI_INTENT_DIMLESS   1011
Dimensionless value - no params - although, as in _ESTIMATE the name of the parameter may be stored in intent_name. 

[DANAJK]

But when you unravel, you need to know the size of each dimension and that is not in the intent code?

I'm not sure what is the point of forcing multi-dimensional data into a single dimension of NIfTI when no-one else will be able to decode the new NIfTI?

What I do in other code is tell my reading programme what parameters I want in dimensions 3+ of my output array. For example, specifying {'slice', 'dynamic'} gives an array that is [nrows ncolumns nslices ndynamics] , or, {'slice','echo'} gives an array that is [nrows ncolumns nslices nechos].

[ashgillman]

But when you unravel, you need to know the size of each dimension and that is not in the intent code?

Yes, that's why it would be encoded it in the intent_name.

I'm not sure what is the point of forcing multi-dimensional data into a single dimension of NIfTI when no-one else will be able to decode the new NIfTI?

I thought the issue was we have only 3 dims to play with in NIfti (5, 6, 7) and 4 dimensions to stuff (channel, phase, reps, sets)?
No matter what is chosen, it will be a unique-to-SIRF encoding.

What I do in other code is tell my reading programme what parameters I want in dimensions 3+ of my output array. For example, specifying {'slice', 'dynamic'} gives an array that is [nrows ncolumns nslices ndynamics] , or, {'slice','echo'} gives an array that is [nrows ncolumns nslices nechos].

Assuming that you'll never need more than 3 dims (does seem quite unlikely to me?) then yes I think this makes more sense to me, to put a key like that into intent_name and have it code what is in dims 5, 6, 7.

[KrisThielemans]

Thanks @ashgillman. Great link to the Nifti doc. Just to repeat

NIFTI-1 reserves dimensions 1,2,3 for space (x,y,z), 4 for time (t), and 5,6,7 for anything else needed.

All solutions will be unsatisfactory. This is true for PET as well. Nifti doesn't store all we need to know (it doesn't even store time frame info). Something to remember is that a SIRF NiftiImageData object could store more meta-info than what is in the nifti header. Output could then be a "meta-data header" + nii. This is for instance done in BIDS, and we have a proposal for STIR with Interfile headers. I guess this could even be used to internally have higher dimensional objects in sirf::NiftiImageData, but limit to 8 dims when writing to file, although I'm sure that will confuse everybody. If someone knows how this is handled in BIDS, maybe that could help.

Given our rate of deciding/implementing this kind of thing, a simple solution is needed though.

I personally have no idea how we currently handle this in MRImageData, but would it be easy enough to let NiftiImageData handle 7D data (of which 3 space and 1 time), where we use dim 5,6,7 according to whatever is needed, and "invent" our own INTENT_CODE for that (e.g. "sirf_phase", "sirf_channel" "sirf_gate" etc). We can then throw an error if we need more than 7.

@ckolbitsch, tagging you as well...

[DANAJK]

But when you unravel, you need to know the size of each dimension and that is not in the intent code?

Yes, that's why it would be encoded it in the intent_name.

Sorry - I mis-understood.

What I do in other code is tell my reading programme what parameters I want in dimensions 3+ of my output array. For example, specifying {'slice', 'dynamic'} gives an array that is [nrows ncolumns nslices ndynamics] , or, {'slice','echo'} gives an array that is [nrows ncolumns nslices nechos].

Assuming that you'll never need more than 3 dims (does seem quite unlikely to me?) then yes I think this makes more sense to me, to put a key like that into intent_name and have it code what is in dims 5, 6, 7.

I don't limit the number of dimensions. This was just an example.