giftiio.write deprecation; suggested alternative misleading

[tomfrankkirk]

Nibabel version 2.5.1

I have a gifti functional image object that I would like to save. When I use top level save function, as in nibabel.save(obj, 'gifti.func.gii'), the files that are produced are not readable by other tools (though nibabel does not give any warnings).

When I use giftiio.write(obj, 'gifti.func.gii'), the produced files interact nicely with other tools, though I get the following deprecation warning DeprecationWarning: giftiio.write function deprecated. Use nibabel.load() instead.

There are two issues: nibabel.save() does not produce useable functional giftis, unlike giftiio.write(), and the suggested alternative of nibabel.load() does not make sense in this context anyway.

[effigies]

This is very strange, because giftiio.write() just passes the image and filename to nibabel.save():

nibabel/nibabel/gifti/giftiio.py

Lines 36 to 85 in a616f7f

	@deprecate_with_version('giftiio.write function deprecated. '
	"Use nibabel.load() instead.",
	'2.1', '4.0')
	def write(image, filename):
	""" Save the current image to a new file
	Parameters
	----------
	image : GiftiImage
	A GiftiImage instance to store
	filename : string
	Filename to store the Gifti file to
	Returns
	-------
	None
	Notes
	-----
	We write all files with utf-8 encoding, and specify this at the top of the
	XML file with the ``encoding`` attribute.
	The Gifti spec suggests using the following suffixes to your
	filename when saving each specific type of data:
	.gii
	Generic GIFTI File
	.coord.gii
	Coordinates
	.func.gii
	Functional
	.label.gii
	Labels
	.rgba.gii
	RGB or RGBA
	.shape.gii
	Shape
	.surf.gii
	Surface
	.tensor.gii
	Tensors
	.time.gii
	Time Series
	.topo.gii
	Topology
	The Gifti file is stored in endian convention of the current machine.
	"""
	from ..loadsave import save
	return save(image, filename)

Could you provide an example that generates this error?

Agreed that it should say save, rather than load...

[tomfrankkirk]

Sure thing - will cook one up tomorrow sometime.

[tomfrankkirk]

Update: i think I have tracked this down to me cutting corners when creating the gifti object. At any rate I can't seem to reproduce this on my linux machine, if I can do so on OS X (where it first came up) then I will put the code here.

[effigies]

Understandable. GIFTI is a pretty complicated format. We recently (#793) improved the API for fetching data from images, and are open to improving the API for creating well-formed images, as well, in case you have any ideas. (cc @htwangtw)

[tomfrankkirk]

I do have one suggestion, but I know its not very good practice... that said, many people do it, because it's an easy way to make things work! When I'm being lazy, I'll load in an existing gifti and simply over-write the data arrays and save it under a new name. Maybe there could a place for a function along the lines of 'create gifti from existing': that uses an existing well-formed gifti as the basis to create a new one, pulling in the appropriate parts of the header. For example, if I've calculated some functional data on a surface, I could use this function along with the surface to create a well-formed func.gii to store the data. What do you think?

[htwangtw]

Thanks @effigies for including me in this thread!
The suggestion from @tomfrankkirk sounds like a practice carried from nifti file creation (copy header, replace data etc.), which I am sure that we all do from time to time. However, this is not as easy to implement in nifti because a gifti file is consist of many gifti arrays and each of them has a header.
From my understanding of gifti files, the crucial information to make a valid file are:

1. data
2. datatype (something like NIFTI_TYPE_FLOAT32 )
   The rest of the header information is handle by some default setting to my knowledge.

I cannot access the code for now, but it's something along this line for the timeseries:

new_gii = empty gifti image
nii_datatype = some valid gifti datatype
data = an array of number of vertices x number of timeseries 
numDA = number of data array; in this case number of timeseries 
for i in numDA:
    new_gii.add_gifti_data_array(GiftiDataArray(data[:, i], nii_datatype ))

Maybe we can have an example on how to create image or a function?
Let me know what I can help.

[tomfrankkirk]

I think an example would certainly be helpful. I've done some more digging today and worked out that type-casting is also giving me grief. The following email thread is similar to the prolem I was having: https://mail.python.org/pipermail//neuroimaging/2017-April/001381.html

Basically I was creating GiftiDataArrays from float64 np arrays, then saving them and having trouble reading them with wb_command. When I explcitly cast them to float32 on creation the problem goes away. Granted, I don't think this is a nibabel problem but its very annoying.

One other comment: the doc page for GiftiDataArray makes it look like the NIFTI intent code should be a string along the lines of NIFTI_INTENT_FLOAT32 but the specified datatype is actually int - bit confusing? The function signature also reads the default as intent='NIFTI_INTENT_NONE'. Reproduced below:

datatype int:
NIFTI_TYPE_UINT8 : Unsigned, 8-bit bytes. NIFTI_TYPE_INT32 : Signed, 32-bit integers. NIFTI_TYPE_FLOAT32 : 32-bit single precision floating point.”

[effigies]

When I explcitly cast them to float32 on creation the problem goes away. Granted, I don't think this is a nibabel problem but its very annoying.

This is definitely frustrating. The question we have to ask is whether coercing to the standard is worth potentially surprising people with lost precision. In general nibabel takes a "programmer knows best" approach, and leaves it to downstream libraries/apps to do coercion.

On the other hand, I don't think any of the binary formats accept data types that the standard does not explicitly permit, so it's probably sensible to move to restrict. We currently coerce at serialization time to the specified NIFTI_TYPE_*, although its default value is just the dtype of the input array.

It would be worth taking a survey of how the binary formats behave when an unsupported dtype is passed, to try to achieve some consistency.

One other comment: the doc page for GiftiDataArray makes it look like the NIFTI intent code should be a string along the lines of NIFTI_INTENT_FLOAT32 but the specified datatype is actually int - bit confusing?

Yes, this could probably be written better. What it means is that you can use any of the equivalent representations in the data type codes, and it will be stored internally as the integer representation:

nibabel/nibabel/nifti1.py

Line 96 in 95c934b

(16, 'float32', np.float32, "NIFTI_TYPE_FLOAT32"),

>>> import numpy as np
>>> from nibabel.nifti1 import data_type_codes
>>> data_type_codes.code['NIFTI_TYPE_FLOAT32']
16
>>> data_type_codes.code['float32']
16
>>> data_type_codes.code[np.float32]
16
>>> data_type_codes.code[16]
16

I'm honestly not sure why we use the int representation as opposed to str. I guess comparisons are cheaper?

[htwangtw]

Took me some time to dig out the code. @tomfrankkirk
Here's a minimum example to create some random values on the 32k gifti HCP surface

import nibabel as nb
import numpy as np
# generic HCP 32k surf left hemi
surf_nii = nb.load('Q1-Q6_R440.L.midthickness.32k_fs_LR.surf.gii')  
coords, triangels = surf_nii.agg_data() 
val = np.random.randn(coords.shape[0]).astype(np.float32)
lh_val_gii = nb.gifti.GiftiImage()                                     
lh_val_gii.add_gifti_data_array(nb.gifti.GiftiDataArray(val))        
lh_val_gii.to_filename('test_32k.L.func.gii') 

When loading with Workbench, select the valid strucure lable (CortexLeft)

[tomfrankkirk]

Thanks for the explanation @effigies and example @htwangtw - I didn't know about the agg_data() shortcut to get the points and triangles out of a surface gifti, thats very handy to know! I understand your general point about trusting that the programmer knows best. Sneaky though some of the issues highlighted here are, I guess they are just part of the learning curve and hopefully a few hints and examples in the docs could help people avoid stumbling on them in the future!