Precision Saving uint16 array

[stefan-kuczera-gu]

Reading in a file with uint16 as data type, saving it again and then reading it again results inconsistent data values.
In the example below an integer value of 3 is converted to 2.999999999301508. Converting the data array from get_fdata() from float64 to uint16 solves the problem. Is this working as intended? In my opinion the conversion should not be necessary. This can be an issue when segmentation masks are edited. File attached below.

import numpy as np
import nibabel as nb

filepath='myseg.nii.gz'
data_obj=nb.load(filepath)
print('obj',data_obj.get_data_dtype())
data=data_obj.get_fdata()
print('numpy arr',np.max(data),data.dtype)
new_manualseg = nb.Nifti1Image(data,affine=data_obj.affine, header=data_obj.header)
nb.save(new_manualseg,'test.nii.gz',dtype=np.uint16)
data=np.array(data,dtype=np.uint16)
new_manualseg = nb.Nifti1Image(data,affine=data_obj.affine, header=data_obj.header)
nb.save(new_manualseg,'test2.nii.gz')

filepath='test.nii.gz'
data_obj=nb.load(filepath)
print('obj',data_obj.get_data_dtype())
data=data_obj.get_fdata()
diff=np.count_nonzero(data-np.rint(data))
print('numpy arr',np.max(data),diff,data.dtype)

filepath='test2.nii.gz'
data_obj=nb.load(filepath)
print('obj',data_obj.get_data_dtype())
data=data_obj.get_fdata()
diff=np.count_nonzero(data-np.rint(data))
print('numpy arr',np.max(data),diff,data.dtype)

Output

obj uint16
numpy arr 3.0 float64
obj uint16
numpy arr 2.999999999301508 347 float64
obj uint16
numpy arr 3.0 0 float64

myseg.nii.gz

[effigies]

Generally, if you pass in a float64 data array and save it as uint or int, then nibabel will attempt to rescale the data so the full range is used. In general, it's an expensive check to verify that all values in an array are perfectly representable as the target type, so that's left to the user.

Anyway, get_fdata() is intended for most MRI data, which samples a continuous range, however it's been discretized. Segmentations are a bit different, and make sense to retain as integer. For reading a segmentation, you'll want to access the dataobj directly:

img = nb.load(filepath)
data = np.asanyarray(img.dataobj)  # Returns the raw data array if no scale factors, float64 if scale factors
data = np.asanyarray(img.dataobj, dtype='uint16')  # Returns uint16 array, regardless of on-disk type or scale factors
data = np.uint16(img.dataobj)  # Same as above

In all cases, the data is scaled if scale factors are set, but with the latter two methods, you have control over the output type.

For creating new arrays where the input type might be different, I would be explicit:

new_img = nb.Nifti1Image(data_array.astype(target_type), img.affine, img.header)

If you really want to pass your data in as floats and not scale, you can also manually set the slope and intercept:

new_img = nb.Nifti1Image(data_array, img.affine, img.header)
new_img.header.set_slope_inter(slope=1, inter=0)