#!/usr/bin/env fslpython import shlex import subprocess as sp import numpy as np import nibabel as nib from fsl.data.image import Image import fsl.transform.affine as affine def applywarp(src, ref): cmd = f'applywarp --in={src} --ref={ref} --usesqform --out=out' sp.run(shlex.split(cmd), check=True) return Image('out') def make_image(fname, data, sform, sform_code, qform, qform_code): hdr = nib.Nifti1Header() hdr.set_sform(sform, sform_code) hdr.set_qform(qform, qform_code) img = Image(data, header=hdr) img.save(fname) return img # fsl/fugue!6 # # Make sure that applywarp --usesqform # will correctly use the sform or qform # # FSL affine precedence rules are as follows: # # 1. If sform_code != 0, use sform; else # 2. If qform_code != 0, use qform; else # 3. Use a scaling matrix with pixdims along the # diagonal def test_applywarp_uses_sqform_correctly(): srcdata = np.zeros((20, 20, 20)) refdata = np.zeros((20, 20, 20)) # Data from the two images # aligned via their affines - # src should be resampled so # as to be identical to ref srcdata[ 6: 9, 6: 9, 6: 9] = 1 refdata[16:19, 16:19, 16:19] = 1 eye = np.eye(4) srcaff = affine.scaleOffsetXform(1, [ 5, 5, 5]) refaff = affine.scaleOffsetXform(1, [-5, -5, -5]) # aligned via sform src = make_image('src', srcdata, srcaff, 2, eye, 1) ref = make_image('ref', refdata, refaff, 2, eye, 1) result = applywarp('src', 'ref') assert np.all(np.isclose(result.data, ref.data)) # aligned via sform again src = make_image('src', srcdata, srcaff, 2, eye, 0) ref = make_image('ref', refdata, refaff, 2, eye, 0) result = applywarp('src', 'ref') assert np.all(np.isclose(result.data, ref.data)) # aligned via qform src = make_image('src', srcdata, eye, 0, srcaff, 2) ref = make_image('ref', refdata, eye, 0, refaff, 2) result = applywarp('src', 'ref') assert np.all(np.isclose(result.data, ref.data)) # not aligned - should result in a scaling # matrix being used (an identify matrix here) src = make_image('src', srcdata, eye, 0, srcaff, 0) ref = make_image('ref', refdata, eye, 0, refaff, 0) result = applywarp('src', 'ref') assert not np.all(np.isclose(result.data, ref.data)) assert np.all(np.isclose(result.data, src.data)) if __name__ == '__main__': test_applywarp_uses_sqform_correctly()