From cf340457b51359efae403e92ad655a3630540905 Mon Sep 17 00:00:00 2001 From: Ben Acland Date: Wed, 26 Feb 2014 10:09:03 -0600 Subject: [PATCH 1/6] some imports --- gradunwarp/core/unwarp_resample.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/gradunwarp/core/unwarp_resample.py b/gradunwarp/core/unwarp_resample.py index 6ee7859..bbaa06e 100644 --- a/gradunwarp/core/unwarp_resample.py +++ b/gradunwarp/core/unwarp_resample.py @@ -7,6 +7,7 @@ import numpy as np import sys import pdb +import gc import math import logging from scipy import ndimage @@ -16,6 +17,7 @@ import globals from globals import siemens_max_det import nibabel as nib +import subprocess #np.seterr(all='raise') From b9826cc92a58791940487df28e8bbbffe3846f2b Mon Sep 17 00:00:00 2001 From: Ben Acland Date: Wed, 26 Feb 2014 10:22:27 -0600 Subject: [PATCH 2/6] add optional fileName argument --- gradunwarp/core/gradient_unwarp.py | 2 +- gradunwarp/core/unwarp_resample.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/gradunwarp/core/gradient_unwarp.py b/gradunwarp/core/gradient_unwarp.py index ce94932..948988e 100644 --- a/gradunwarp/core/gradient_unwarp.py +++ b/gradunwarp/core/gradient_unwarp.py @@ -96,7 +96,7 @@ def run(self): self.vol, self.m_rcs2ras = utils.get_vol_affine(self.args.infile) self.unwarper = Unwarper(self.vol, self.m_rcs2ras, self.args.vendor, - self.coeffs) + self.coeffs, getattr(self.args, 'infile', None)) if hasattr(self.args, 'fovmin') and self.args.fovmin: self.unwarper.fovmin = float(self.args.fovmin) if hasattr(self.args, 'fovmax') and self.args.fovmax: diff --git a/gradunwarp/core/unwarp_resample.py b/gradunwarp/core/unwarp_resample.py index bbaa06e..70477bc 100644 --- a/gradunwarp/core/unwarp_resample.py +++ b/gradunwarp/core/unwarp_resample.py @@ -27,13 +27,14 @@ class Unwarper(object): ''' ''' - def __init__(self, vol, m_rcs2ras, vendor, coeffs): + def __init__(self, vol, m_rcs2ras, vendor, coeffs, fileName=None): ''' ''' self.vol = vol self.m_rcs2ras = m_rcs2ras self.vendor = vendor self.coeffs = coeffs + self.name=fileName self.warp = False self.nojac = False self.m_rcs2lai = None @@ -95,7 +96,6 @@ def eval_spharm_grid(self, vendor, coeffs): return CV(_dv.x, _dv.y, _dv.z), CV(gr, gc, gs), g_xyz2rcs - def run(self): ''' ''' From 57a8d724bf820ca06cd27bf12485cfe8bc474d66 Mon Sep 17 00:00:00 2001 From: Ben Acland Date: Wed, 26 Feb 2014 10:50:06 -0600 Subject: [PATCH 3/6] always convert to float32 if out dtype is float64 --- gradunwarp/core/unwarp_resample.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/gradunwarp/core/unwarp_resample.py b/gradunwarp/core/unwarp_resample.py index 70477bc..25f20e6 100644 --- a/gradunwarp/core/unwarp_resample.py +++ b/gradunwarp/core/unwarp_resample.py @@ -181,10 +181,12 @@ def run(self): def write(self, outfile): log.info('Writing output to ' + outfile) + # if out datatype is float64 make it float32 + if self.out.dtype == np.float64: + self.out = self.out.astype(np.float32) if outfile.endswith('.nii') or outfile.endswith('.nii.gz'): img = nib.Nifti1Image(self.out, self.m_rcs2ras) if outfile.endswith('.mgh') or outfile.endswith('.mgz'): - self.out = self.out.astype(np.float32) img = nib.MGHImage(self.out, self.m_rcs2ras) nib.save(img, outfile) From 6840664a2979ff44fe75fad7974d5bed29b7e627 Mon Sep 17 00:00:00 2001 From: Ben Acland Date: Wed, 26 Feb 2014 10:51:21 -0600 Subject: [PATCH 4/6] get out of the way, print --- gradunwarp/core/unwarp_resample.py | 1 - 1 file changed, 1 deletion(-) diff --git a/gradunwarp/core/unwarp_resample.py b/gradunwarp/core/unwarp_resample.py index 25f20e6..8772a92 100644 --- a/gradunwarp/core/unwarp_resample.py +++ b/gradunwarp/core/unwarp_resample.py @@ -172,7 +172,6 @@ def run(self): dvz[..., slice] = moddv.z ''' - print # Evaluate spherical harmonics on a smaller grid dv, grcs, g_xyz2rcs = self.eval_spharm_grid(self.vendor, self.coeffs) # do the nonlinear unwarp From fbcbe2c410684c605478dde7d57bc6e98c3010b3 Mon Sep 17 00:00:00 2001 From: Ben Acland Date: Wed, 26 Feb 2014 11:09:54 -0600 Subject: [PATCH 5/6] use non_linear_unwarp_siemens if possible, and requested added in the --hcp argument to avoid overloading filename. It's possible that the two non_linear_unwarp methods could be merged in the future, but I didn't write either, and don't feel like fiddling with them. So here. --- gradunwarp/core/gradient_unwarp.py | 4 +- gradunwarp/core/unwarp_resample.py | 190 +++++++++++++++++++++++++++-- 2 files changed, 183 insertions(+), 11 deletions(-) mode change 100644 => 100755 gradunwarp/core/unwarp_resample.py diff --git a/gradunwarp/core/gradient_unwarp.py b/gradunwarp/core/gradient_unwarp.py index 948988e..e07cb6b 100644 --- a/gradunwarp/core/gradient_unwarp.py +++ b/gradunwarp/core/gradient_unwarp.py @@ -47,6 +47,8 @@ def argument_parse_gradunwarp(): help='number of grid points in each direction') p.add_argument('--interp_order', dest='order', help='the order of interpolation(1..4) where 1 is linear - default') + p.add_argument('--hcp', dest='hcp', action='store_true', default=False, + help='if true, and infile is provided, and vendor is siemens, try to use non_linear_unwarp_siemens.') p.add_argument('--verbose', action='store_true', default=False) @@ -96,7 +98,7 @@ def run(self): self.vol, self.m_rcs2ras = utils.get_vol_affine(self.args.infile) self.unwarper = Unwarper(self.vol, self.m_rcs2ras, self.args.vendor, - self.coeffs, getattr(self.args, 'infile', None)) + self.coeffs, self.args.infile, self.args.hcp) if hasattr(self.args, 'fovmin') and self.args.fovmin: self.unwarper.fovmin = float(self.args.fovmin) if hasattr(self.args, 'fovmax') and self.args.fovmax: diff --git a/gradunwarp/core/unwarp_resample.py b/gradunwarp/core/unwarp_resample.py old mode 100644 new mode 100755 index 8772a92..d849ade --- a/gradunwarp/core/unwarp_resample.py +++ b/gradunwarp/core/unwarp_resample.py @@ -27,7 +27,7 @@ class Unwarper(object): ''' ''' - def __init__(self, vol, m_rcs2ras, vendor, coeffs, fileName=None): + def __init__(self, vol, m_rcs2ras, vendor, coeffs, fileName=None, try_hcp=False): ''' ''' self.vol = vol @@ -38,6 +38,7 @@ def __init__(self, vol, m_rcs2ras, vendor, coeffs, fileName=None): self.warp = False self.nojac = False self.m_rcs2lai = None + self.try_hcp = try_hcp # grid is uninit by default self.fovmin = None @@ -105,6 +106,9 @@ def run(self): if self.warp: self.polarity = -1. + # use the HCP siemens unwarp method if possible + use_hcp_siemens = (self.try_hcp and self.vendor == 'siemens' and self.name != None) + # transform RAS-coordinates into LAI-coordinates m_ras2lai = np.array([[-1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], @@ -113,16 +117,18 @@ def run(self): m_rcs2lai = np.dot(m_ras2lai, self.m_rcs2ras) # indices of image volume - nr, nc, ns = self.vol.shape[:3] - vc3, vr3, vs3 = utils.meshgrid(np.arange(nr), np.arange(nc), np.arange(ns), dtype=np.float32) - vrcs = CV(x=vr3, y=vc3, z=vs3) - vxyz = utils.transform_coordinates(vrcs, m_rcs2lai) + if not use_hcp_siemens: + nr, nc, ns = self.vol.shape[:3] + vc3, vr3, vs3 = utils.meshgrid(np.arange(nr), np.arange(nc), np.arange(ns), dtype=np.float32) + vrcs = CV(x=vr3, y=vc3, z=vs3) + vxyz = utils.transform_coordinates(vrcs, m_rcs2lai) # account for half-voxel shift in R and C directions - halfvox = np.zeros((4, 4)) - halfvox[0, 3] = m_rcs2lai[0, 0] / 2.0 - halfvox[1, 3] = m_rcs2lai[1, 1] / 2.0 - m_rcs2lai = m_rcs2lai + halfvox + if not use_hcp_siemens: + halfvox = np.zeros((4, 4)) + halfvox[0, 3] = m_rcs2lai[0, 0] / 2.0 + halfvox[1, 3] = m_rcs2lai[1, 1] / 2.0 + m_rcs2lai = m_rcs2lai + halfvox # extract rotational and scaling parts of the transformation matrix # ignore the translation part @@ -175,7 +181,11 @@ def run(self): # Evaluate spherical harmonics on a smaller grid dv, grcs, g_xyz2rcs = self.eval_spharm_grid(self.vendor, self.coeffs) # do the nonlinear unwarp - self.out, self.vjacmult_lps = self.non_linear_unwarp(vxyz, grcs, dv, dxyz, + if use_hcp_siemens: + self.out, self.vjacout = self.non_linear_unwarp_siemens(self.vol.shape, dv, dxyz, + m_rcs2lai, m_rcs2lai, g_xyz2rcs) + else: + self.out, self.vjacmult_lps = self.non_linear_unwarp(vxyz, grcs, dv, dxyz, m_rcs2lai, g_xyz2rcs) def write(self, outfile): @@ -302,6 +312,166 @@ def non_linear_unwarp(self, vxyz, grcs, dv, dxyz, m_rcs2lai, g_xyz2rcs): if self.vendor == 'ge': pass # for now + def non_linear_unwarp_siemens(self, volshape, dv, dxyz, m_rcs2lai, m_rcs2lai, g_xyz2rcs): + ''' Performs the crux of the unwarping. + It's agnostic to Siemens or GE and uses more functions to + do the processing separately. + + Needs self.vendor, self.coeffs, self.warp, self.nojac to be set + + Parameters + ---------- + vxyz : CoordsVector (namedtuple) contains np.array + has 3 elements x,y and z each representing the grid coordinates + dxyz : CoordsVector (namedtuple) + differential coords vector + + Returns + ------- + TODO still vague what to return + vwxyz : CoordsVector (namedtuple) contains np.array + x,y and z coordinates of the unwarped coordinates + vjacmult_lps : np.array + the jacobian multiplier (determinant) + ''' + log.info('Evaluating the jacobian multiplier') + nr, nc, ns = self.vol.shape[:3] + if not self.nojac: + jim2 = np.zeros((nr, nc), dtype=np.float32) + vjacdet_lpsw = np.zeros((nr, nc), dtype=np.float32) + if dxyz == 0: + vjacdet_lps = 1 + else: + vjacdet_lps = eval_siemens_jacobian_mult(dv, dxyz) + + # essentially pre-allocating everything + out = np.zeros((nr, nc, ns), dtype=np.float32) + fullWarp = np.zeros((nr, nc, ns, 3), dtype=np.float32) + + vjacout = np.zeros((nr, nc, ns), dtype=np.float32) + im2 = np.zeros((nr, nc), dtype=np.float32) + dvx = np.zeros((nr, nc), dtype=np.float32) + dvy = np.zeros((nr, nc), dtype=np.float32) + dvz = np.zeros((nr, nc), dtype=np.float32) + im_ = np.zeros((nr, nc), dtype=np.float32) + # init jacobian temp image + vc, vr = utils.meshgrid(np.arange(nc), np.arange(nr)) + + log.info('Unwarping slice by slice') + # for every slice + for s in xrange(ns): + # pretty print + sys.stdout.flush() + if (s+1) % 10 == 0: + print s+1, + else: + print '.', + + # hopefully, free memory + gc.collect() + # init to 0 + dvx.fill(0.) + dvy.fill(0.) + dvz.fill(0.) + im_.fill(0.) + + vs = np.ones(vr.shape) * s + vrcs = CV(vr, vc, vs) + vxyz = utils.transform_coordinates(vrcs, m_rcs2lai) + vrcsg = utils.transform_coordinates(vxyz, g_xyz2rcs) + ndimage.interpolation.map_coordinates(dv.x, + vrcsg, + output=dvx, + order=self.order) + ndimage.interpolation.map_coordinates(dv.y, + vrcsg, + output=dvy, + order=self.order) + ndimage.interpolation.map_coordinates(dv.z, + vrcsg, + output=dvz, + order=self.order) + # new locations of the image voxels in XYZ ( LAI ) coords + + #dvx.fill(0.) + #dvy.fill(0.) + #dvz.fill(0.) + + vxyzw = CV(x=vxyz.x + self.polarity * dvx, + y=vxyz.y + self.polarity * dvy, + z=vxyz.z + self.polarity * dvz) + + # convert the locations got into RCS indices + vrcsw = utils.transform_coordinates(vxyzw, + np.linalg.inv(m_rcs2lai)) + # map the internal voxel coordinates to FSL scaled mm coordinates + pixdim1=float((subprocess.Popen(['fslval', self.name,'pixdim1'], stdout=subprocess.PIPE).communicate()[0]).strip()) + pixdim2=float((subprocess.Popen(['fslval', self.name,'pixdim2'], stdout=subprocess.PIPE).communicate()[0]).strip()) + pixdim3=float((subprocess.Popen(['fslval', self.name,'pixdim3'], stdout=subprocess.PIPE).communicate()[0]).strip()) + dim1=float((subprocess.Popen(['fslval', self.name,'dim1'], stdout=subprocess.PIPE).communicate()[0]).strip()) + outputOrient=subprocess.Popen(['fslorient', self.name], stdout=subprocess.PIPE).communicate()[0] + if outputOrient.strip() == 'NEUROLOGICAL': + # if neurological then flip x coordinate (both here in premat and later in postmat) + m_vox2fsl = np.array([[-1.0*pixdim1, 0.0, 0.0, pixdim1*(dim1-1)], + [0.0, pixdim2, 0.0, 0.0], + [0.0, 0.0, pixdim3, 0.0], + [0.0, 0.0, 0.0, 1.0]], dtype=np.float) + else: + m_vox2fsl = np.array([[pixdim1, 0.0, 0.0, 0.0], + [0.0, pixdim2, 0.0, 0.0], + [0.0, 0.0, pixdim3, 0.0], + [0.0, 0.0, 0.0, 1.0]], dtype=np.float) + + vfsl = utils.transform_coordinates(vrcsw, m_vox2fsl) + + + #im_ = utils.interp3(self.vol, vrcsw.x, vrcsw.y, vrcsw.z) + ndimage.interpolation.map_coordinates(self.vol, + vrcsw, + output=im_, + order=self.order) + # find NaN voxels, and set them to 0 + im_[np.where(np.isnan(im_))] = 0. + im_[np.where(np.isinf(im_))] = 0. + im2[vr, vc] = im_ + + #img = nib.Nifti1Image(dvx,np.eye(4)) + #nib.save(img,"x"+str(s).zfill(3)+".nii.gz") + #img = nib.Nifti1Image(dvy,np.eye(4)) + #nib.save(img,"y"+str(s).zfill(3)+".nii.gz") + #img = nib.Nifti1Image(dvz,np.eye(4)) + #nib.save(img,"z"+str(s).zfill(3)+".nii.gz") + + # Multiply the intensity with the Jacobian det, if needed + if not self.nojac: + vjacdet_lpsw.fill(0.) + jim2.fill(0.) + # if polarity is negative, the jacobian is also inversed + if self.polarity == -1: + vjacdet_lps = 1. / vjacdet_lps + + ndimage.interpolation.map_coordinates(vjacdet_lps, + vrcsg, + output=vjacdet_lpsw, + order=self.order) + vjacdet_lpsw[np.where(np.isnan(vjacdet_lpsw))] = 0. + vjacdet_lpsw[np.where(np.isinf(vjacdet_lpsw))] = 0. + jim2[vr, vc] = vjacdet_lpsw + im2 = im2 * jim2 + vjacout[..., s] = jim2 + + fullWarp[...,s,0]=vfsl.x + fullWarp[...,s,1]=vfsl.y + fullWarp[...,s,2]=vfsl.z + out[..., s] = im2 + + print + + img=nib.Nifti1Image(fullWarp,self.m_rcs2ras) + nib.save(img,"fullWarp_abs.nii.gz") + # return image and the jacobian + return out, vjacout + def eval_siemens_jacobian_mult(F, dxyz): ''' From cb5cc5d83dc553271abbeccf7cd4a3359126ca19 Mon Sep 17 00:00:00 2001 From: Ben Acland Date: Wed, 26 Feb 2014 11:21:57 -0600 Subject: [PATCH 6/6] remove redundant arg So, in the original hcp modifications, there was this m_rcs2lai_nohalf variable... but the line that applied halfvox to m_rcs2lai was also commented out. Now the halfvox stuff only happens when use_hcp_siemens is false, and there never was any difference between m_rcs2lai and m_rcs2lai_nohalf. So I renamed all _nohalf to remove the suffix, but didn't notice that the two were being redundantly passed into use_hcp_siemens. --- gradunwarp/core/unwarp_resample.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/gradunwarp/core/unwarp_resample.py b/gradunwarp/core/unwarp_resample.py index d849ade..e16c350 100755 --- a/gradunwarp/core/unwarp_resample.py +++ b/gradunwarp/core/unwarp_resample.py @@ -183,7 +183,7 @@ def run(self): # do the nonlinear unwarp if use_hcp_siemens: self.out, self.vjacout = self.non_linear_unwarp_siemens(self.vol.shape, dv, dxyz, - m_rcs2lai, m_rcs2lai, g_xyz2rcs) + m_rcs2lai, g_xyz2rcs) else: self.out, self.vjacmult_lps = self.non_linear_unwarp(vxyz, grcs, dv, dxyz, m_rcs2lai, g_xyz2rcs) @@ -312,7 +312,7 @@ def non_linear_unwarp(self, vxyz, grcs, dv, dxyz, m_rcs2lai, g_xyz2rcs): if self.vendor == 'ge': pass # for now - def non_linear_unwarp_siemens(self, volshape, dv, dxyz, m_rcs2lai, m_rcs2lai, g_xyz2rcs): + def non_linear_unwarp_siemens(self, volshape, dv, dxyz, m_rcs2lai, g_xyz2rcs): ''' Performs the crux of the unwarping. It's agnostic to Siemens or GE and uses more functions to do the processing separately.