DerivativesDataSink inconsistent byte order

[liningpan]

What happened?

DerivativesDataSink could produce invalid/corrupted NIfTI file, where the endianness of the header is inconsistent with the data volume.

fmriprep use DerivativesDataSink to coerce averaged boldref into float32 data type. If the input NIfTI file is in big endian and the volume require data type conversion (i.e. in i16), the call to unsafe write will generate an output nifti file with a big endian header and a little endian (machine native) data volume, which will obviously corrupt downstream steps. (PS: antsAI will seg fault if none of the start locations worked, which made debugging this even less obvious).

niworkflows/niworkflows/interfaces/bids.py

Line 1152 in 9e18356

data_dtype = self.inputs.data_dtype or self._default_dtypes[self.inputs.suffix]

niworkflows/niworkflows/interfaces/bids.py

Lines 1210 to 1223 in 9e18356

	if new_data is new_header is None:
	_copy_any(orig_file, str(out_file))
	else:
	orig_img = nb.load(orig_file)
	if new_data is None:
	set_consumables(new_header, orig_img.dataobj)
	new_data = orig_img.dataobj.get_unscaled()
	else:
	# Without this, we would be writing nans
	# This is our punishment for hacking around nibabel defaults
	new_header.set_slope_inter(slope=1.0, inter=0.0)
	unsafe_write_nifti_header_and_data(
	fname=out_file, header=new_header, data=new_data
	)

niworkflows/niworkflows/utils/images.py

Line 52 in 9e18356

def unsafe_write_nifti_header_and_data(fname, header, data):

What command did you use?

I ran into this issue through fmriprep. This is a intrinsic operation in niworkflows

What version of the software are you running?

fmriprep 25.2.3 sigularity (built from docker)

How are you running this software?

Singularity

Is your data BIDS valid?

Yes

Are you reusing any previously computed results?

No

Please copy and paste any relevant log output.

Additional information / screenshots

No response

[oesteban]

cc/ @effigies What do you think?

[effigies]

I don't think I've ever seen a big-endian NIfTI, so that would explain why we never saw this. @liningpan Could you create a minimal reproducible example that we could use as a test?

[liningpan]

Here is a minimal reproducible example. This is derived from intermediate results of https://openneuro.org/datasets/ds000258, many but not all scans in this dataset are in big-endian format. I don't think I've seen a big-endian NIfTI file in the wild before either.

from niworkflows.interfaces.bids import DerivativesDataSink


def main():
    ds = DerivativesDataSink(
        in_file="sub-21281_task-rest_echo-1_bold_average.nii",
        source_file="sub-21281/func/sub-21281_task-rest_echo-1_bold.nii",
        base_directory=".",
        desc="hmc",
        suffix="boldref",
        compress=False,
        dismiss_entities=["echo"],
    )
    ds.run()

    print("input")
    with open("sub-21281_task-rest_echo-1_bold_average.nii", "rb") as f:
        data = f.read(352 + 16)
        print_hex(data)
    """input
    0000 00 00 01 5c 00 00 00 00
    0008 00 00 00 00 00 00 00 00
    0010 00 00 00 00 00 00 00 00
    0018 00 00 00 00 00 00 00 00
    0020 00 00 40 00 00 00 72 00
    0028 00 03 00 40 00 40 00 1e
    0030 00 01 00 01 00 01 00 01
    0038 00 00 00 00 00 00 00 00
    0040 00 00 00 00 00 00 00 04
    0048 00 10 00 00 3f 80 00 00
    0050 40 70 00 00 40 70 00 00
    0058 40 8c cc cd 3f 80 00 00
    0060 3f 80 00 00 3f 80 00 00
    0068 3f 80 00 00 43 b0 00 00
    0070 3c f1 a4 4e 44 71 a4 4e <- scl_slope 0.029497292 (BE) scl_intercept 966.56726
    0078 00 1d 00 0a 44 c2 e0 00
    0080 00 00 00 00 00 00 00 00
    0088 00 00 00 00 00 00 00 00
    0090 00 00 00 00 00 00 00 00
    0098 00 00 00 00 00 00 00 00
    00a0 00 00 00 00 00 00 00 00
    00a8 00 00 00 00 00 00 00 00
    00b0 00 00 00 00 00 00 00 00
    00b8 00 00 00 00 00 00 00 00
    00c0 00 00 00 00 00 00 00 00
    00c8 00 00 00 00 00 00 00 00
    00d0 00 00 00 00 00 00 00 00
    00d8 00 00 00 00 00 00 00 00
    00e0 00 00 00 00 00 00 00 00
    00e8 00 00 00 00 00 00 00 00
    00f0 00 00 00 00 00 00 00 00
    00f8 00 00 00 00 00 01 00 01
    0100 3c a3 4f 1d be 81 62 ce
    0108 3f 77 a2 f0 42 f7 5f 18
    0110 43 15 19 05 41 cd 46 f7
    0118 c0 6f cd 63 bd 85 f1 c8
    0120 3e 25 1a 18 42 f7 5f 18
    0128 bc 26 86 2b c0 51 57 01
    0130 c0 09 b7 b1 43 15 19 05
    0138 3e 1b 7d 26 bf ea 99 a0
    0140 40 75 68 e8 41 cd 46 f7
    0148 00 00 00 00 00 00 00 00
    0150 00 00 00 00 00 00 00 00
    0158 6e 2b 31 00 00 00 00 00
    0160 80 00 81 29 81 43 81 3e <- -32768 -32471 -32445 -32450 (Big-endian i2) 
                                    apply nifti scale 0, 8.76069, 9.52762, 9.38013
    """

    print("output")
    with open(
        "niworkflows/sub-21281/func/sub-21281_task-rest_desc-hmc_boldref.nii", "rb"
    ) as f:
        data = f.read(352 + 16)
        print_hex(data)

    """output
    0000 00 00 01 5c 00 00 00 00 <- header is still big-endian
    0008 00 00 00 00 00 00 00 00
    0010 00 00 00 00 00 00 00 00
    0018 00 00 00 00 00 00 00 00
    0020 00 00 40 00 00 00 72 00
    0028 00 03 00 40 00 40 00 1e
    0030 00 01 00 01 00 01 00 01
    0038 00 00 00 00 00 00 00 00
    0040 00 00 00 00 00 00 00 10
    0048 00 20 00 00 3f 80 00 00
    0050 40 70 00 00 40 70 00 00
    0058 40 8c cc cd 3f 80 00 00
    0060 3f 80 00 00 3f 80 00 00 <- scl slope and intercept are 1
    0068 3f 80 00 00 43 b0 00 00
    0070 3f 80 00 00 00 00 00 00
    0078 00 1d 00 02 44 c2 e0 00
    0080 00 00 00 00 00 00 00 00
    0088 00 00 00 00 00 00 00 00
    0090 00 00 00 00 00 00 00 00
    0098 00 00 00 00 00 00 00 00
    00a0 00 00 00 00 00 00 00 00
    00a8 00 00 00 00 00 00 00 00
    00b0 00 00 00 00 00 00 00 00
    00b8 00 00 00 00 00 00 00 00
    00c0 00 00 00 00 00 00 00 00
    00c8 00 00 00 00 00 00 00 00
    00d0 00 00 00 00 00 00 00 00
    00d8 00 00 00 00 00 00 00 00
    00e0 00 00 00 00 00 00 00 00
    00e8 00 00 00 00 00 00 00 00
    00f0 00 00 00 00 00 00 00 00
    00f8 00 00 00 00 00 01 00 01
    0100 3c a3 4f 1c be 81 62 ce
    0108 3f 77 a2 f0 42 f7 5f 18
    0110 43 15 19 05 41 cd 46 f7
    0118 c0 6f cd 63 bd 85 f1 c8
    0120 3e 25 1a 18 42 f7 5f 18
    0128 bc 26 86 2b c0 51 57 01
    0130 c0 09 b7 b1 43 15 19 05
    0138 3e 1b 7d 26 bf ea 99 a0
    0140 40 75 68 e8 41 cd 46 f7
    0148 00 00 00 00 00 00 00 00
    0150 00 00 00 00 00 00 00 00
    0158 6e 2b 31 00 00 00 00 00
    0160 00 00 00 00 cf 2b 0c 41 <- 0, 8.760695 (in little endian f32)
    0168 27 71 18 41 0c 15 16 41 <- 9.527625, 9.380138
    """


def print_hex(data):
    for i in range(0, len(data), 8):
        chunk = data[i : i + 8]
        print("{:04x} {}".format(i, " ".join(f"{b:02x}" for b in chunk)))


if __name__ == "__main__":
    main()

repo.zip

[effigies]

Thank you!