From 6c84ea8e4884f509ab2e8be5dd76895190fab05e Mon Sep 17 00:00:00 2001 From: mathiasg Date: Fri, 18 Oct 2019 12:27:38 -0400 Subject: [PATCH 01/15] enh: lta file read/write from nipy/nibabel#565 --- nitransforms/io.py | 232 ++++++++++++++++++++++++++++++++ nitransforms/linear.py | 31 +++-- nitransforms/tests/data/inv.lta | 29 ++++ 3 files changed, 284 insertions(+), 8 deletions(-) create mode 100644 nitransforms/io.py create mode 100644 nitransforms/tests/data/inv.lta diff --git a/nitransforms/io.py b/nitransforms/io.py new file mode 100644 index 00000000..9cdf6e45 --- /dev/null +++ b/nitransforms/io.py @@ -0,0 +1,232 @@ +import numpy as np +from nibabel.wrapstruct import LabeledWrapStruct +from nibabel.volumeutils import Recoder + +transform_codes = Recoder(( + (0, 'LINEAR_VOX_TO_VOX'), + (1, 'LINEAR_RAS_TO_RAS'), + (2, 'LINEAR_PHYSVOX_TO_PHYSVOX'), + (14, 'REGISTER_DAT'), + (21, 'LINEAR_COR_TO_COR')), + fields=('code', 'label')) + + +class StringBasedStruct(LabeledWrapStruct): + def __init__(self, + binaryblock=None, + endianness=None, + check=True): + if binaryblock is not None and getattr(binaryblock, 'dtype', + None) == self.dtype: + self._structarr = binaryblock.copy() + return + super(StringBasedStruct, self).__init__(binaryblock, endianness, check) + + def __array__(self): + return self._structarr + + +class VolumeGeometry(StringBasedStruct): + template_dtype = np.dtype([ + ('valid', 'i4'), # Valid values: 0, 1 + ('volume', 'i4', (3, 1)), # width, height, depth + ('voxelsize', 'f4', (3, 1)), # xsize, ysize, zsize + ('xras', 'f4', (3, 1)), # x_r, x_a, x_s + ('yras', 'f4', (3, 1)), # y_r, y_a, y_s + ('zras', 'f4', (3, 1)), # z_r, z_a, z_s + ('cras', 'f4', (3, 1)), # c_r, c_a, c_s + ('filename', 'U1024')]) # Not conformant (may be >1024 bytes) + dtype = template_dtype + + def as_affine(self): + affine = np.eye(4) + sa = self.structarr + A = np.hstack((sa['xras'], sa['yras'], sa['zras'])) * sa['voxelsize'] + b = sa['cras'] - A.dot(sa['volume']) / 2 + affine[:3, :3] = A + affine[:3, [3]] = b + return affine + + def to_string(self): + sa = self.structarr + lines = [ + 'valid = {} # volume info {:s}valid'.format( + sa['valid'], '' if sa['valid'] else 'in'), + 'filename = {}'.format(sa['filename']), + 'volume = {:d} {:d} {:d}'.format(*sa['volume'].flatten()), + 'voxelsize = {:.15e} {:.15e} {:.15e}'.format( + *sa['voxelsize'].flatten()), + 'xras = {:.15e} {:.15e} {:.15e}'.format(*sa['xras'].flatten()), + 'yras = {:.15e} {:.15e} {:.15e}'.format(*sa['yras'].flatten()), + 'zras = {:.15e} {:.15e} {:.15e}'.format(*sa['zras'].flatten()), + 'cras = {:.15e} {:.15e} {:.15e}'.format(*sa['cras'].flatten()), + ] + return '\n'.join(lines) + + @classmethod + def from_image(klass, img): + volgeom = klass() + sa = volgeom.structarr + sa['valid'] = 1 + sa['volume'][:, 0] = img.shape[:3] # Assumes xyzt-ordered image + sa['voxelsize'][:, 0] = img.header.get_zooms()[:3] + A = img.affine[:3, :3] + b = img.affine[:3, [3]] + cols = A * (1 / sa['voxelsize']) + sa['xras'] = cols[:, [0]] + sa['yras'] = cols[:, [1]] + sa['zras'] = cols[:, [2]] + sa['cras'] = b + A.dot(sa['volume']) / 2 + try: + sa['filename'] = img.file_map['image'].filename + except Exception: + pass + + return volgeom + + @classmethod + def from_string(klass, string): + volgeom = klass() + sa = volgeom.structarr + lines = string.splitlines() + for key in ('valid', 'filename', 'volume', 'voxelsize', + 'xras', 'yras', 'zras', 'cras'): + label, valstring = lines.pop(0).split(' = ') + assert label.strip() == key + + val = np.genfromtxt([valstring.encode()], + dtype=klass.dtype[key]) + sa[key] = val.reshape(sa[key].shape) if val.size else '' + + return volgeom + + +class LinearTransform(StringBasedStruct): + template_dtype = np.dtype([ + ('mean', 'f4', (3, 1)), # x0, y0, z0 + ('sigma', 'f4'), + ('m_L', 'f4', (4, 4)), + ('m_dL', 'f4', (4, 4)), + ('m_last_dL', 'f4', (4, 4)), + ('src', VolumeGeometry), + ('dst', VolumeGeometry), + ('label', 'i4')]) + dtype = template_dtype + + def __getitem__(self, idx): + val = super(LinearTransform, self).__getitem__(idx) + if idx in ('src', 'dst'): + val = VolumeGeometry(val) + return val + + def to_string(self): + sa = self.structarr + lines = [ + 'mean = {:6.4f} {:6.4f} {:6.4f}'.format( + *sa['mean'].flatten()), + 'sigma = {:6.4f}'.format(float(sa['sigma'])), + '1 4 4', + ('{:18.15e} ' * 4).format(*sa['m_L'][0]), + ('{:18.15e} ' * 4).format(*sa['m_L'][1]), + ('{:18.15e} ' * 4).format(*sa['m_L'][2]), + ('{:18.15e} ' * 4).format(*sa['m_L'][3]), + 'src volume info', + self['src'].to_string(), + 'dst volume info', + self['dst'].to_string(), + ] + return '\n'.join(lines) + + @classmethod + def from_string(klass, string): + lt = klass() + sa = lt.structarr + lines = string.splitlines() + for key in ('mean', 'sigma'): + label, valstring = lines.pop(0).split(' = ') + assert label.strip() == key + + val = np.genfromtxt([valstring.encode()], + dtype=klass.dtype[key]) + sa[key] = val.reshape(sa[key].shape) + assert lines.pop(0) == '1 4 4' + val = np.genfromtxt([valstring.encode() for valstring in lines[:4]], + dtype='f4') + sa['m_L'] = val + lines = lines[4:] + assert lines.pop(0) == 'src volume info' + sa['src'] = np.asanyarray(VolumeGeometry.from_string('\n'.join(lines[:8]))) + lines = lines[8:] + assert lines.pop(0) == 'dst volume info' + sa['dst'] = np.asanyarray(VolumeGeometry.from_string('\n'.join(lines))) + return lt + + +class LinearTransformArray(StringBasedStruct): + template_dtype = np.dtype([ + ('type', 'i4'), + ('nxforms', 'i4'), + ('subject', 'U1024'), + ('fscale', 'f4')]) + dtype = template_dtype + _xforms = None + + def __init__(self, + binaryblock=None, + endianness=None, + check=True): + super(LinearTransformArray, self).__init__(binaryblock, endianness, check) + self._xforms = [LinearTransform() + for _ in range(self.structarr['nxforms'])] + + def __getitem__(self, idx): + if idx == 'xforms': + return self._xforms + if idx == 'nxforms': + return len(self._xforms) + return super(LinearTransformArray, self).__getitem__(idx) + + def to_string(self): + code = int(self['type']) + header = [ + 'type = {} # {}'.format(code, transform_codes.label[code]), + 'nxforms = {}'.format(self['nxforms'])] + xforms = [xfm.to_string() for xfm in self._xforms] + footer = [ + 'subject {}'.format(self['subject']), + 'fscale {:.6f}'.format(float(self['fscale']))] + return '\n'.join(header + xforms + footer) + + @classmethod + def from_string(klass, string): + lta = klass() + sa = lta.structarr + lines = string.splitlines() + for key in ('type', 'nxforms'): + label, valstring = lines.pop(0).split(' = ') + assert label.strip() == key + + val = np.genfromtxt([valstring.encode()], + dtype=klass.dtype[key]) + sa[key] = val.reshape(sa[key].shape) if val.size else '' + for _ in range(sa['nxforms']): + lta._xforms.append( + LinearTransform.from_string('\n'.join(lines[:25]))) + lines = lines[25:] + for key in ('subject', 'fscale'): + # Optional keys + if not lines[0].startswith(key): + continue + label, valstring = lines.pop(0).split(' ') + assert label.strip() == key + + val = np.genfromtxt([valstring.encode()], + dtype=klass.dtype[key]) + sa[key] = val.reshape(sa[key].shape) if val.size else '' + + assert len(lta._xforms) == sa['nxforms'] + return lta + + @classmethod + def from_fileobj(klass, fileobj, check=True): + return klass.from_string(fileobj.read()) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index b61578ea..676400fb 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -11,11 +11,13 @@ import numpy as np from scipy import ndimage as ndi from pathlib import Path +import warnings from nibabel.loadsave import load as loadimg from nibabel.affines import from_matvec, voxel_sizes, obliquity from .base import TransformBase from .patched import shape_zoom_affine +from .io import LinearTransformArray, VolumeGeometry LPS = np.diag([-1, -1, 1, 1]) @@ -126,8 +128,7 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, try: reference = self.reference except ValueError: - print('Warning: no reference space defined, using moving as reference', - file=sys.stderr) + warnings.warn('No reference space defined, using moving as reference') reference = moving nvols = 1 @@ -270,13 +271,13 @@ def to_filename(self, filename, fmt='X5', moving=None): 3dvolreg matrices (DICOM-to-DICOM, row-by-row):""", fmt='%g') return filename - if fmt.lower() == 'fsl': - if not moving: - moving = self.reference - - if isinstance(moving, str): - moving = loadimg(moving) + # for FSL / FS information + if not moving: + moving = self.reference + elif isinstance(moving, str): + moving = loadimg(moving) + if fmt.lower() == 'fsl': # Adjust for reference image offset and orientation refswp, refspc = _fsl_aff_adapt(self.reference) pre = self.reference.affine.dot( @@ -298,6 +299,14 @@ def to_filename(self, filename, fmt='X5', moving=None): else: np.savetxt(filename, mat[0], delimiter=' ', fmt='%g') return filename + elif fmt.lower() in ('fs', 'lta'): + # linear tranform info + src = VolumeGeometry.from_image(moving) + dst = VolumeGeometry.from_image(self.reference) + + + + return super(Affine, self).to_filename(filename, fmt=fmt) @@ -326,6 +335,12 @@ def load(filename, fmt='X5', reference=None): # elif fmt.lower() == 'afni': # parameters = LPS.dot(self.matrix.dot(LPS)) # parameters = parameters[:3, :].reshape(-1).tolist() + elif fmt.lower() in ('fs', 'lta'): + with open(filename) as ltafile: + lta = LinearTransformArray.from_fileobj(ltafile) + # TODO: multiple transforms? + assert lta['nxforms'] == 1 + matrix = lta._xforms[0]['m_L'] elif fmt.lower() in ('x5', 'bids'): raise NotImplementedError else: diff --git a/nitransforms/tests/data/inv.lta b/nitransforms/tests/data/inv.lta new file mode 100644 index 00000000..68093d48 --- /dev/null +++ b/nitransforms/tests/data/inv.lta @@ -0,0 +1,29 @@ +type = 1 # LINEAR_RAS_TO_RAS +nxforms = 1 +mean = 0.0000 0.0000 0.0000 +sigma = 1.0000 +1 4 4 +9.719529747962952e-01 -2.037279307842255e-02 -8.194014430046082e-03 -1.919340729713440e+00 +1.478777173906565e-02 8.941915035247803e-01 -4.111929237842560e-01 2.351728630065918e+01 +3.128148615360260e-02 4.410418868064880e-01 7.873729467391968e-01 1.280669403076172e+01 +0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +src volume info +valid = 1 # volume info valid +filename = /opt/freesurfer/average/mni305.cor.mgz +volume = 256 256 256 +voxelsize = 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00 +xras = -1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 +yras = 0.000000000000000e+00 0.000000000000000e+00 -1.000000000000000e+00 +zras = 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 +cras = 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 +dst volume info +valid = 1 # volume info valid +filename = /home/jovyan/data/derivatives/mindboggle/freesurfer_subjects/sub-voice969/mri/orig_nu.mgz +volume = 256 256 256 +voxelsize = 1.000000000000000e+00 1.000000000000000e+00 1.000000000000000e+00 +xras = -1.000000000000000e+00 1.396983861923218e-09 0.000000000000000e+00 +yras = -9.313225746154785e-10 0.000000000000000e+00 -9.999998807907104e-01 +zras = 9.313225746154785e-10 1.000000000000000e+00 -2.980232238769531e-08 +cras = -1.583099365234375e-02 3.479890441894531e+01 -6.033630371093750e+00 +subject fsaverage +fscale 0.100000 From e821b80e3a64bc070b8453e87cb8b179c358b44e Mon Sep 17 00:00:00 2001 From: mathiasg Date: Fri, 18 Oct 2019 12:32:02 -0400 Subject: [PATCH 02/15] tst: lta file parsing --- nitransforms/tests/test_io.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) create mode 100644 nitransforms/tests/test_io.py diff --git a/nitransforms/tests/test_io.py b/nitransforms/tests/test_io.py new file mode 100644 index 00000000..fc037900 --- /dev/null +++ b/nitransforms/tests/test_io.py @@ -0,0 +1,23 @@ +import os + +import numpy as np + +from ..io import LinearTransformArray as LTA + + +def test_LinearTransformArray_input(tmpdir, data_path): + lta = LTA() + assert lta['nxforms'] == 0 + assert len(lta['xforms']) == 0 + + test_lta = os.path.join(data_path, 'inv.lta') + with open(test_lta) as fp: + lta = LTA.from_fileobj(fp) + + assert lta.get('type') == 1 + assert len(lta['xforms']) == lta['nxforms'] == 1 + xform = lta['xforms'][0] + + assert np.allclose( + xform['m_L'], np.genfromtxt(test_lta, skip_header=5, skip_footer=20) + ) From 515fb0a5cf5e1b13caf0d297be43778e274d7079 Mon Sep 17 00:00:00 2001 From: mathiasg Date: Fri, 18 Oct 2019 18:14:45 -0400 Subject: [PATCH 03/15] enh: write lta --- nitransforms/io.py | 9 +++++++-- nitransforms/linear.py | 23 +++++++++++++++-------- 2 files changed, 22 insertions(+), 10 deletions(-) diff --git a/nitransforms/io.py b/nitransforms/io.py index 9cdf6e45..929d2673 100644 --- a/nitransforms/io.py +++ b/nitransforms/io.py @@ -1,6 +1,7 @@ import numpy as np from nibabel.wrapstruct import LabeledWrapStruct from nibabel.volumeutils import Recoder +from nibabel.affines import voxel_sizes transform_codes = Recoder(( (0, 'LINEAR_VOX_TO_VOX'), @@ -69,10 +70,10 @@ def from_image(klass, img): sa = volgeom.structarr sa['valid'] = 1 sa['volume'][:, 0] = img.shape[:3] # Assumes xyzt-ordered image - sa['voxelsize'][:, 0] = img.header.get_zooms()[:3] + sa['voxelsize'][:, 0] = voxel_sizes(img.affine)[:3] A = img.affine[:3, :3] b = img.affine[:3, [3]] - cols = A * (1 / sa['voxelsize']) + cols = A / sa['voxelsize'] sa['xras'] = cols[:, [0]] sa['yras'] = cols[:, [1]] sa['zras'] = cols[:, [2]] @@ -208,6 +209,10 @@ def from_string(klass, string): val = np.genfromtxt([valstring.encode()], dtype=klass.dtype[key]) + if val not in (0, 1): + raise NotImplementedError( + "LTA type {0} is not currently supported".format(transform_codes.label[val]) + ) sa[key] = val.reshape(sa[key].shape) if val.size else '' for _ in range(sa['nxforms']): lta._xforms.append( diff --git a/nitransforms/linear.py b/nitransforms/linear.py index 676400fb..4753ebe5 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -17,7 +17,7 @@ from nibabel.affines import from_matvec, voxel_sizes, obliquity from .base import TransformBase from .patched import shape_zoom_affine -from .io import LinearTransformArray, VolumeGeometry +from .io import LinearTransformArray, LinearTransform, VolumeGeometry LPS = np.diag([-1, -1, 1, 1]) @@ -300,12 +300,20 @@ def to_filename(self, filename, fmt='X5', moving=None): np.savetxt(filename, mat[0], delimiter=' ', fmt='%g') return filename elif fmt.lower() in ('fs', 'lta'): - # linear tranform info - src = VolumeGeometry.from_image(moving) - dst = VolumeGeometry.from_image(self.reference) - - + # xform info + lt = LinearTransform() + lt['sigma'] = 1. + lt['m_L'] = self.matrix + lt['src'] = VolumeGeometry.from_image(moving) + lt['dst'] = VolumeGeometry.from_image(self.reference) + # to make LTA file format + lta = LinearTransformArray() + lta['type'] = 1 # RAS2RAS + lta['xforms'] = [lt] + with open(filename, 'w') as f: + f.write(lta.to_string()) + return filename return super(Affine, self).to_filename(filename, fmt=fmt) @@ -338,8 +346,7 @@ def load(filename, fmt='X5', reference=None): elif fmt.lower() in ('fs', 'lta'): with open(filename) as ltafile: lta = LinearTransformArray.from_fileobj(ltafile) - # TODO: multiple transforms? - assert lta['nxforms'] == 1 + assert lta['nxforms'] == 1 # ever have multiple transforms? matrix = lta._xforms[0]['m_L'] elif fmt.lower() in ('x5', 'bids'): raise NotImplementedError From 9557221b73f779a79b9ad239995b879dae182647 Mon Sep 17 00:00:00 2001 From: mathiasg Date: Sat, 19 Oct 2019 13:21:35 -0400 Subject: [PATCH 04/15] enh: interal vox2ras conversion --- nitransforms/io.py | 52 +++++++++++++++++++++++++++++++----------- nitransforms/linear.py | 9 ++++---- 2 files changed, 44 insertions(+), 17 deletions(-) diff --git a/nitransforms/io.py b/nitransforms/io.py index 929d2673..69f45044 100644 --- a/nitransforms/io.py +++ b/nitransforms/io.py @@ -209,25 +209,22 @@ def from_string(klass, string): val = np.genfromtxt([valstring.encode()], dtype=klass.dtype[key]) - if val not in (0, 1): - raise NotImplementedError( - "LTA type {0} is not currently supported".format(transform_codes.label[val]) - ) sa[key] = val.reshape(sa[key].shape) if val.size else '' for _ in range(sa['nxforms']): lta._xforms.append( LinearTransform.from_string('\n'.join(lines[:25]))) lines = lines[25:] - for key in ('subject', 'fscale'): - # Optional keys - if not lines[0].startswith(key): - continue - label, valstring = lines.pop(0).split(' ') - assert label.strip() == key + if lines: + for key in ('subject', 'fscale'): + # Optional keys + if not lines[0].startswith(key): + continue + label, valstring = lines.pop(0).split(' ') + assert label.strip() == key - val = np.genfromtxt([valstring.encode()], - dtype=klass.dtype[key]) - sa[key] = val.reshape(sa[key].shape) if val.size else '' + val = np.genfromtxt([valstring.encode()], + dtype=klass.dtype[key]) + sa[key] = val.reshape(sa[key].shape) if val.size else '' assert len(lta._xforms) == sa['nxforms'] return lta @@ -235,3 +232,32 @@ def from_string(klass, string): @classmethod def from_fileobj(klass, fileobj, check=True): return klass.from_string(fileobj.read()) + + def as_type(self, target): + """ + Convert the internal transformation matrix to a different type inplace + + Parameters + ---------- + target : str, int + Tranformation type + """ + assert self['nxforms'] == 1, "Cannot convert multiple transformations" + xform = self['xforms'][0] + src = xform['src'] + dst = xform['dst'] + current = self['type'] + if isinstance(target, str): + target = transform_codes.code[target] + + # VOX2VOX -> RAS2RAS + if current == 0 and target == 1: + M = np.linalg.inv(src.as_affine()).dot(xform['m_L']).dot(dst.as_affine()) + xform['m_L'] = M + else: + raise NotImplementedError( + "Converting {0} to {1} is not yet available".format( + transform_codes.label[current], + transform_codes.label[target] + ) + ) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index 4753ebe5..ffc8220f 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -151,8 +151,7 @@ def resample(self, moving, order=3, mode='constant', cval=0.0, prefilter=True, singlemat = np.linalg.inv(movaff).dot(self._matrix[0].dot(reference.affine)) if singlemat is not None and nvols > nmats: - print('Warning: resampling a 4D volume with a single affine matrix', - file=sys.stderr) + warnings.warn('Resampling a 4D volume with a single affine matrix') # Compose an index to index affine matrix moved = [] @@ -274,7 +273,7 @@ def to_filename(self, filename, fmt='X5', moving=None): # for FSL / FS information if not moving: moving = self.reference - elif isinstance(moving, str): + if isinstance(moving, str): moving = loadimg(moving) if fmt.lower() == 'fsl': @@ -347,7 +346,9 @@ def load(filename, fmt='X5', reference=None): with open(filename) as ltafile: lta = LinearTransformArray.from_fileobj(ltafile) assert lta['nxforms'] == 1 # ever have multiple transforms? - matrix = lta._xforms[0]['m_L'] + if lta['type'] != 1: + lta.as_type(1) + matrix = lta['xforms'][0]['m_L'] elif fmt.lower() in ('x5', 'bids'): raise NotImplementedError else: From 7d1362ba3add05e19254dc431805c2259db49e7e Mon Sep 17 00:00:00 2001 From: Mathias Goncalves Date: Sat, 19 Oct 2019 15:58:37 -0400 Subject: [PATCH 05/15] Update nitransforms/io.py Co-Authored-By: Chris Markiewicz --- nitransforms/io.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nitransforms/io.py b/nitransforms/io.py index 69f45044..9af304d8 100644 --- a/nitransforms/io.py +++ b/nitransforms/io.py @@ -233,7 +233,7 @@ def from_string(klass, string): def from_fileobj(klass, fileobj, check=True): return klass.from_string(fileobj.read()) - def as_type(self, target): + def set_type(self, target): """ Convert the internal transformation matrix to a different type inplace From dc6c6b19410bf906284593d2a874b885e5886f8b Mon Sep 17 00:00:00 2001 From: mathiasg Date: Mon, 21 Oct 2019 14:46:31 -0400 Subject: [PATCH 06/15] fix: change LTA type --- nitransforms/io.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/nitransforms/io.py b/nitransforms/io.py index 9af304d8..77bf85d3 100644 --- a/nitransforms/io.py +++ b/nitransforms/io.py @@ -253,7 +253,6 @@ def set_type(self, target): # VOX2VOX -> RAS2RAS if current == 0 and target == 1: M = np.linalg.inv(src.as_affine()).dot(xform['m_L']).dot(dst.as_affine()) - xform['m_L'] = M else: raise NotImplementedError( "Converting {0} to {1} is not yet available".format( @@ -261,3 +260,5 @@ def set_type(self, target): transform_codes.label[target] ) ) + xform['m_L'] = M + self['type'] = target From 8e23723ddaa0fd67457d1e33307b736672627854 Mon Sep 17 00:00:00 2001 From: mathiasg Date: Mon, 21 Oct 2019 14:48:05 -0400 Subject: [PATCH 07/15] fix: alter lta type, only use lta extension --- nitransforms/linear.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index ffc8220f..739029bd 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -298,7 +298,7 @@ def to_filename(self, filename, fmt='X5', moving=None): else: np.savetxt(filename, mat[0], delimiter=' ', fmt='%g') return filename - elif fmt.lower() in ('fs', 'lta'): + elif fmt.lower() == 'lta': # xform info lt = LinearTransform() lt['sigma'] = 1. @@ -342,12 +342,13 @@ def load(filename, fmt='X5', reference=None): # elif fmt.lower() == 'afni': # parameters = LPS.dot(self.matrix.dot(LPS)) # parameters = parameters[:3, :].reshape(-1).tolist() - elif fmt.lower() in ('fs', 'lta'): + elif fmt.lower() == 'lta': with open(filename) as ltafile: lta = LinearTransformArray.from_fileobj(ltafile) - assert lta['nxforms'] == 1 # ever have multiple transforms? + if lta['nxforms'] > 1: + raise NotImplementedError("Multiple transforms are not yet supported.") if lta['type'] != 1: - lta.as_type(1) + lta.set_type(1) matrix = lta['xforms'][0]['m_L'] elif fmt.lower() in ('x5', 'bids'): raise NotImplementedError From d0a05765fa68fc4bedd77347ae11347241a23ddd Mon Sep 17 00:00:00 2001 From: mathiasg Date: Mon, 21 Oct 2019 14:48:38 -0400 Subject: [PATCH 08/15] enh: LTA + components tests --- nitransforms/tests/test_io.py | 37 +++++++++++++++++++++++++++++++++-- 1 file changed, 35 insertions(+), 2 deletions(-) diff --git a/nitransforms/tests/test_io.py b/nitransforms/tests/test_io.py index fc037900..ebf1ce5f 100644 --- a/nitransforms/tests/test_io.py +++ b/nitransforms/tests/test_io.py @@ -2,10 +2,35 @@ import numpy as np -from ..io import LinearTransformArray as LTA +from ..io import ( + VolumeGeometry as VG, + LinearTransform as LT, + LinearTransformArray as LTA, +) -def test_LinearTransformArray_input(tmpdir, data_path): +def test_VolumeGeometry(tmpdir, get_data): + vg = VG() + assert vg['valid'] == 0 + + img = get_data['RAS'] + vg = VG.from_image(img) + assert vg['valid'] == 1 + assert np.all(vg['voxelsize'] == img.header.get_zooms()[:3]) + assert np.all(vg.as_affine() == img.affine) + + assert len(vg.to_string().split('\n')) == 8 + + +def test_LinearTransform(tmpdir, get_data): + lt = LT() + assert lt['m_L'].shape == (4, 4) + assert np.all(lt['m_L'] == 0) + for vol in ('src', 'dst'): + assert lt[vol]['valid'] == 0 + + +def test_LinearTransformArray(tmpdir, data_path): lta = LTA() assert lta['nxforms'] == 0 assert len(lta['xforms']) == 0 @@ -21,3 +46,11 @@ def test_LinearTransformArray_input(tmpdir, data_path): assert np.allclose( xform['m_L'], np.genfromtxt(test_lta, skip_header=5, skip_footer=20) ) + + outlta = (tmpdir / 'out.lta').strpath + with open(outlta, 'w') as fp: + fp.write(lta.to_string()) + + with open(outlta) as fp: + lta2 = LTA.from_fileobj(fp) + np.allclose(lta['xforms'][0]['m_L'], lta2['xforms'][0]['m_L']) From 4b18c65d1c5a46d72c142f164bbdcdf638168cae Mon Sep 17 00:00:00 2001 From: mathiasg Date: Mon, 21 Oct 2019 16:28:57 -0400 Subject: [PATCH 09/15] fix: patch warp struct --- nitransforms/io.py | 5 +++-- nitransforms/patched.py | 14 ++++++++++++++ 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/nitransforms/io.py b/nitransforms/io.py index 77bf85d3..d45d1d67 100644 --- a/nitransforms/io.py +++ b/nitransforms/io.py @@ -1,8 +1,9 @@ import numpy as np -from nibabel.wrapstruct import LabeledWrapStruct from nibabel.volumeutils import Recoder from nibabel.affines import voxel_sizes +from .patched import LabeledWrapStruct + transform_codes = Recoder(( (0, 'LINEAR_VOX_TO_VOX'), (1, 'LINEAR_RAS_TO_RAS'), @@ -150,7 +151,7 @@ def from_string(klass, string): val = np.genfromtxt([valstring.encode()], dtype=klass.dtype[key]) sa[key] = val.reshape(sa[key].shape) - assert lines.pop(0) == '1 4 4' + assert lines.pop(0) == '1 4 4' # xforms, shape + 1, shape + 1 val = np.genfromtxt([valstring.encode() for valstring in lines[:4]], dtype='f4') sa['m_L'] = val diff --git a/nitransforms/patched.py b/nitransforms/patched.py index 6e0b7ecb..1bfbea9c 100644 --- a/nitransforms/patched.py +++ b/nitransforms/patched.py @@ -1,4 +1,5 @@ import numpy as np +from nibabel.wrapstruct import LabeledWrapStruct as LWS def shape_zoom_affine(shape, zooms, x_flip=True, y_flip=False): @@ -63,3 +64,16 @@ def shape_zoom_affine(shape, zooms, x_flip=True, y_flip=False): aff[:3, :3] = np.diag(zooms) aff[:3, -1] = -origin * zooms return aff + + +class LabeledWrapStruct(LWS): + def __setitem__(self, item, value): + ''' Set values in structured data + Examples + -------- + >>> wstr = WrapStruct() + >>> wstr['integer'] = 3 + >>> wstr['integer'] + array(3, dtype=int16) + ''' + self._structarr[item] = np.asanyarray(value) From bf2e90b28787d7c1202060cdce4301f925e966b2 Mon Sep 17 00:00:00 2001 From: mathiasg Date: Mon, 21 Oct 2019 16:29:53 -0400 Subject: [PATCH 10/15] fix: LTA file write --- nitransforms/linear.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index 739029bd..d342206d 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -17,7 +17,7 @@ from nibabel.affines import from_matvec, voxel_sizes, obliquity from .base import TransformBase from .patched import shape_zoom_affine -from .io import LinearTransformArray, LinearTransform, VolumeGeometry +from . import io LPS = np.diag([-1, -1, 1, 1]) @@ -300,15 +300,15 @@ def to_filename(self, filename, fmt='X5', moving=None): return filename elif fmt.lower() == 'lta': # xform info - lt = LinearTransform() + lt = io.LinearTransform() lt['sigma'] = 1. lt['m_L'] = self.matrix - lt['src'] = VolumeGeometry.from_image(moving) - lt['dst'] = VolumeGeometry.from_image(self.reference) + lt['src'] = io.VolumeGeometry.from_image(moving) + lt['dst'] = io.VolumeGeometry.from_image(self.reference) # to make LTA file format - lta = LinearTransformArray() + lta = io.LinearTransformArray() lta['type'] = 1 # RAS2RAS - lta['xforms'] = [lt] + lta['xforms'].append(lt) with open(filename, 'w') as f: f.write(lta.to_string()) @@ -344,7 +344,7 @@ def load(filename, fmt='X5', reference=None): # parameters = parameters[:3, :].reshape(-1).tolist() elif fmt.lower() == 'lta': with open(filename) as ltafile: - lta = LinearTransformArray.from_fileobj(ltafile) + lta = io.LinearTransformArray.from_fileobj(ltafile) if lta['nxforms'] > 1: raise NotImplementedError("Multiple transforms are not yet supported.") if lta['type'] != 1: From a0788131ed051473bd0c524713c995141d892f38 Mon Sep 17 00:00:00 2001 From: mathiasg Date: Mon, 21 Oct 2019 16:57:14 -0400 Subject: [PATCH 11/15] enh: linear lta read/write testing --- nitransforms/tests/data/affine-LAS.fs.lta | 29 +++++++++++++++++++++++ nitransforms/tests/data/affine-LPS.fs.lta | 29 +++++++++++++++++++++++ nitransforms/tests/data/affine-RAS.fs.lta | 29 +++++++++++++++++++++++ nitransforms/tests/test_transform.py | 12 +++++++--- 4 files changed, 96 insertions(+), 3 deletions(-) create mode 100644 nitransforms/tests/data/affine-LAS.fs.lta create mode 100644 nitransforms/tests/data/affine-LPS.fs.lta create mode 100644 nitransforms/tests/data/affine-RAS.fs.lta diff --git a/nitransforms/tests/data/affine-LAS.fs.lta b/nitransforms/tests/data/affine-LAS.fs.lta new file mode 100644 index 00000000..7fc69a8a --- /dev/null +++ b/nitransforms/tests/data/affine-LAS.fs.lta @@ -0,0 +1,29 @@ +type = 1 # LINEAR_RAS_TO_RAS +nxforms = 1 +mean = 0.0000 0.0000 0.0000 +sigma = 1.0000 +1 4 4 +9.999989867210388e-01 -9.999993490055203e-04 9.999998146668077e-04 4.000000000000000e+00 +1.404936308972538e-03 6.216088533401489e-01 -7.833265066146851e-01 2.000000000000000e+00 +1.617172238184139e-04 7.833271622657776e-01 6.216096282005310e-01 -1.000000000000000e+00 +0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +src volume info +valid = 1 # volume info valid +filename = +volume = 57 67 56 +voxelsize = 2.750000000000000e+00 2.750000000000000e+00 2.750000000000000e+00 +xras = -1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 +yras = 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 +zras = 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +cras = -2.375000000000000e+00 1.125000000000000e+00 -1.400000000000000e+01 +dst volume info +valid = 1 # volume info valid +filename = +volume = 57 67 56 +voxelsize = 2.750000000000000e+00 2.750000000000000e+00 2.750000000000000e+00 +xras = -1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 +yras = 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 +zras = 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +cras = -2.375000000000000e+00 1.125000000000000e+00 -1.400000000000000e+01 +subject +fscale 0.000000 \ No newline at end of file diff --git a/nitransforms/tests/data/affine-LPS.fs.lta b/nitransforms/tests/data/affine-LPS.fs.lta new file mode 100644 index 00000000..3a9d7c1c --- /dev/null +++ b/nitransforms/tests/data/affine-LPS.fs.lta @@ -0,0 +1,29 @@ +type = 1 # LINEAR_RAS_TO_RAS +nxforms = 1 +mean = 0.0000 0.0000 0.0000 +sigma = 1.0000 +1 4 4 +9.999989867210388e-01 -9.999993490055203e-04 9.999998146668077e-04 4.000000000000000e+00 +1.404936308972538e-03 6.216088533401489e-01 -7.833265066146851e-01 2.000000000000000e+00 +1.617172238184139e-04 7.833271622657776e-01 6.216096282005310e-01 -1.000000000000000e+00 +0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +src volume info +valid = 1 # volume info valid +filename = +volume = 57 67 56 +voxelsize = 2.750000000000000e+00 2.750000000000000e+00 2.750000000000000e+00 +xras = -1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 +yras = 0.000000000000000e+00 -1.000000000000000e+00 0.000000000000000e+00 +zras = 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +cras = -2.375000000000000e+00 -1.625000000000000e+00 -1.400000000000000e+01 +dst volume info +valid = 1 # volume info valid +filename = +volume = 57 67 56 +voxelsize = 2.750000000000000e+00 2.750000000000000e+00 2.750000000000000e+00 +xras = -1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 +yras = 0.000000000000000e+00 -1.000000000000000e+00 0.000000000000000e+00 +zras = 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +cras = -2.375000000000000e+00 -1.625000000000000e+00 -1.400000000000000e+01 +subject +fscale 0.000000 \ No newline at end of file diff --git a/nitransforms/tests/data/affine-RAS.fs.lta b/nitransforms/tests/data/affine-RAS.fs.lta new file mode 100644 index 00000000..eef71d72 --- /dev/null +++ b/nitransforms/tests/data/affine-RAS.fs.lta @@ -0,0 +1,29 @@ +type = 1 # LINEAR_RAS_TO_RAS +nxforms = 1 +mean = 0.0000 0.0000 0.0000 +sigma = 1.0000 +1 4 4 +9.999989867210388e-01 -9.999993490055203e-04 9.999998146668077e-04 4.000000000000000e+00 +1.404936308972538e-03 6.216088533401489e-01 -7.833265066146851e-01 2.000000000000000e+00 +1.617172238184139e-04 7.833271622657776e-01 6.216096282005310e-01 -1.000000000000000e+00 +0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +src volume info +valid = 1 # volume info valid +filename = +volume = 57 67 56 +voxelsize = 2.750000000000000e+00 2.750000000000000e+00 2.750000000000000e+00 +xras = 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 +yras = 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 +zras = 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +cras = 3.750000000000000e-01 1.125000000000000e+00 -1.400000000000000e+01 +dst volume info +valid = 1 # volume info valid +filename = +volume = 57 67 56 +voxelsize = 2.750000000000000e+00 2.750000000000000e+00 2.750000000000000e+00 +xras = 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 +yras = 0.000000000000000e+00 1.000000000000000e+00 0.000000000000000e+00 +zras = 0.000000000000000e+00 0.000000000000000e+00 1.000000000000000e+00 +cras = 3.750000000000000e-01 1.125000000000000e+00 -1.400000000000000e+01 +subject +fscale 0.000000 \ No newline at end of file diff --git a/nitransforms/tests/test_transform.py b/nitransforms/tests/test_transform.py index 23c8db7b..72662cdd 100644 --- a/nitransforms/tests/test_transform.py +++ b/nitransforms/tests/test_transform.py @@ -27,7 +27,7 @@ 'afni': """\ 3dAllineate -base {reference} -input {moving} \ -prefix resampled.nii.gz -1Dmatrix_apply {transform} -final NN\ -""".format, +""".format,- } @@ -35,7 +35,7 @@ @pytest.mark.parametrize('image_orientation', [ 'RAS', 'LAS', 'LPS', # 'oblique', ]) -@pytest.mark.parametrize('sw_tool', ['itk', 'fsl', 'afni']) +@pytest.mark.parametrize('sw_tool', ['itk', 'fsl', 'afni', 'fs']) def test_linear_load(tmpdir, data_path, get_data, image_orientation, sw_tool): """Check implementation of loading affines from formats.""" tmpdir.chdir() @@ -51,6 +51,8 @@ def test_linear_load(tmpdir, data_path, get_data, image_orientation, sw_tool): ext = '' if sw_tool == 'itk': ext = '.tfm' + elif sw_tool == 'fs': + ext = '.lta' fname = 'affine-%s.%s%s' % (image_orientation, sw_tool, ext) xfm_fname = os.path.join(data_path, fname) @@ -74,6 +76,8 @@ def test_linear_load(tmpdir, data_path, get_data, image_orientation, sw_tool): loaded = nbl.load(xfm_fname, fmt=fmt, reference='img.nii.gz') elif sw_tool == 'itk': loaded = nbl.load(xfm_fname, fmt=fmt) + elif sw_tool == 'fs': + loaded = nbl.load(xfm_fname, fmt=fmt) assert loaded == xfm @@ -81,7 +85,7 @@ def test_linear_load(tmpdir, data_path, get_data, image_orientation, sw_tool): @pytest.mark.parametrize('image_orientation', [ 'RAS', 'LAS', 'LPS', # 'oblique', ]) -@pytest.mark.parametrize('sw_tool', ['itk', 'fsl', 'afni']) +@pytest.mark.parametrize('sw_tool', ['itk', 'fsl', 'afni', 'fs']) def test_linear_save(data_path, get_data, image_orientation, sw_tool): """Check implementation of exporting affines to formats.""" img = get_data[image_orientation] @@ -93,6 +97,8 @@ def test_linear_save(data_path, get_data, image_orientation, sw_tool): ext = '' if sw_tool == 'itk': ext = '.tfm' + elif sw_tool == 'fs': + ext = '.lta' with InTemporaryDirectory(): xfm_fname1 = 'M.%s%s' % (sw_tool, ext) From 45be03b700bd528567b0862b3a0990d57c359211 Mon Sep 17 00:00:00 2001 From: mathiasg Date: Mon, 21 Oct 2019 17:00:45 -0400 Subject: [PATCH 12/15] fix: syntax --- nitransforms/tests/test_transform.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/nitransforms/tests/test_transform.py b/nitransforms/tests/test_transform.py index 72662cdd..c7a03b27 100644 --- a/nitransforms/tests/test_transform.py +++ b/nitransforms/tests/test_transform.py @@ -27,7 +27,7 @@ 'afni': """\ 3dAllineate -base {reference} -input {moving} \ -prefix resampled.nii.gz -1Dmatrix_apply {transform} -final NN\ -""".format,- +""".format, } From 3ee09f5f0a70c6c94b602e79b9c092562776f385 Mon Sep 17 00:00:00 2001 From: mathiasg Date: Mon, 21 Oct 2019 18:04:02 -0400 Subject: [PATCH 13/15] fix: fmt type --- nitransforms/linear.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/nitransforms/linear.py b/nitransforms/linear.py index d342206d..9ffd2d0d 100644 --- a/nitransforms/linear.py +++ b/nitransforms/linear.py @@ -298,7 +298,7 @@ def to_filename(self, filename, fmt='X5', moving=None): else: np.savetxt(filename, mat[0], delimiter=' ', fmt='%g') return filename - elif fmt.lower() == 'lta': + elif fmt.lower() == 'fs': # xform info lt = io.LinearTransform() lt['sigma'] = 1. @@ -342,7 +342,7 @@ def load(filename, fmt='X5', reference=None): # elif fmt.lower() == 'afni': # parameters = LPS.dot(self.matrix.dot(LPS)) # parameters = parameters[:3, :].reshape(-1).tolist() - elif fmt.lower() == 'lta': + elif fmt.lower() == 'fs': with open(filename) as ltafile: lta = io.LinearTransformArray.from_fileobj(ltafile) if lta['nxforms'] > 1: From 77fdbc991ec28dd54c648be8e5df520eea4f2de3 Mon Sep 17 00:00:00 2001 From: mathiasg Date: Mon, 21 Oct 2019 18:04:18 -0400 Subject: [PATCH 14/15] enh: support lta loading --- nitransforms/tests/utils.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/nitransforms/tests/utils.py b/nitransforms/tests/utils.py index 3bea4357..cc01a3ec 100644 --- a/nitransforms/tests/utils.py +++ b/nitransforms/tests/utils.py @@ -11,9 +11,15 @@ def assert_affines_by_filename(affine1, affine2): affine2 = Path(affine2) assert affine1.suffix == affine2.suffix, 'affines of different type' - if affine1.suffix.endswith('.tfm'): # An ITK transform - xfm1 = nbl.load(str(affine1), fmt='itk') - xfm2 = nbl.load(str(affine2), fmt='itk') + ext_to_fmt = { + '.tfm': 'itk', # An ITK transform + '.lta': 'fs', # FreeSurfer LTA + } + + ext = affine1.suffix[-4:] + if ext in ext_to_fmt: + xfm1 = nbl.load(str(affine1), fmt=ext_to_fmt[ext]) + xfm2 = nbl.load(str(affine2), fmt=ext_to_fmt[ext]) assert xfm1 == xfm2 else: xfm1 = np.loadtxt(str(affine1)) From 1aee2c4dc81c8b220c6c28eb87bb3e968b898fa1 Mon Sep 17 00:00:00 2001 From: mathiasg Date: Mon, 21 Oct 2019 18:05:51 -0400 Subject: [PATCH 15/15] tst: remove nibabel doctest --- nitransforms/patched.py | 8 -------- 1 file changed, 8 deletions(-) diff --git a/nitransforms/patched.py b/nitransforms/patched.py index 1bfbea9c..ec47a667 100644 --- a/nitransforms/patched.py +++ b/nitransforms/patched.py @@ -68,12 +68,4 @@ def shape_zoom_affine(shape, zooms, x_flip=True, y_flip=False): class LabeledWrapStruct(LWS): def __setitem__(self, item, value): - ''' Set values in structured data - Examples - -------- - >>> wstr = WrapStruct() - >>> wstr['integer'] = 3 - >>> wstr['integer'] - array(3, dtype=int16) - ''' self._structarr[item] = np.asanyarray(value)