BUG: Complex data cast to float in get_fdata()

[amirshamaei]

Hi all, I was wondering how nibabel deals with complex numbers. I tried to store complex data with a nifti library in java, then import it using nibabel. the original data type is complex 64 which got stored in the first dimension(i.e 'x') (e.g imagine data is a time series 10 points and no spatial information; so the storage size is [2,1,1,10]). when I load it in nibabel, it can identify the data type(complex64); but it uses float64 to read the array and I don't have any clue how it reads the imaginary part. I appreciate if anyone gives me some info about the complex number in nifti

[effigies]

This might be a problem where you're using img.get_fdata(), which will coerce to float64 by default, losing the imaginary component. You should be able to use img.get_fdata(dtype=np.complex64) or np.asanyarray(img.dataobj) to get the data as it's encoded on disk.

[amirshamaei]

@effigies thanks. it helped. this issue has been solved. *surprisingly when I am using complex128 the results are accurate as complex64.

[matthew-brett]

Oops - yes - we should probably return complex here - don't you think?

[effigies]

Yeah, I think detecting the complex case is worth doing. I think it's not uncommon to see get_fdata(dtype=np.float32) as people try to avoid a float64 cast; what makes sense to do in that case? TypeError or "you probably know what you're doing"?

[amirshamaei]

sorry for closing the topic :) It would be really nice to detect complex. I think spec2nii API uses nibabel as a backend for accessing nifti files; and spectroscopy imaging data usually has complex numbers.

[matthew-brett]

@effigies - yes I think dtype=np.float32 is a know-what-you're-doing case. At least, if they read their own code, they shouldn't be suprised they got a float32!

[effigies]

Currently the default dtype=np.float64. Going to need to think through a consistent logic to ensure that the input is respected while providing a way to say "give me the narrowest dtype that fits the data" or "give me the narrowest dtype that fits the data and the provided dtype".

For the first case, ideally we'd get:

on-disk dtype	output dtype
int	float32 or float64, depending on values
float32	float32
float64	float64
complex64	complex64

For the latter case:

dtype arg	on-disk dtype	output dtype
float32	int	float32
float64	int	float64
float32	float64	float64
float32	complex64	complex64
float32	complex128	complex128
float64	complex64	complex128

I think the fundamental idea is captured by numpy.result_type, but the result should be unsurprising. Possibly we need a new argument to determine the type resolution mode.

[matthew-brett]

Oh - I was thinking only that, if the image-implied dtype was complex, we'd return complex128 by default, instead of float64, and we'd respect the explicit dtype whatver that was.