Added coil combination using espirit

coil sensitivies can be saved to nifti and loaded back in
enabling them to be calculated on one set of data and used on anouther

concept.py

@@ -14,26 +14,28 @@ import datetime, json, nibabel as nib
 import numpy as np
 import finufft  #flatron insitute nufft
 
+
+
 #local  
-import siemens
+import siemens, espirit
 
 class Recon:
-  """
-  initiate paramters for reconstruction of concept
-  fname : path to twix file
-  optional:
-    alpha=1, scale the FOV 
-             set larger than 1 will reconstruct a larger FOV, the Image matrix  
-             remains the same (i.e. lower resolution)
-    was_alpha_used_in_measurment=False, set to True if sequence .ini file 
-             applied alpha to the chosen FOV (ensure correct FOV size)
-    Nppr=64, default is 64, number RO points pere ring
-    ImgRescale=1 rescale the image matrix (i.e. resolution)
-             set to 2 to double reconstructed matrix size, maintaining the FOV
-  """
-  
-  def __init__(self,fname,alpha=1.,Nppr=64,ImgRescale=1,was_alpha_used_in_measurment=False,is_density_weighted=True,filter_kspace=True):
-    self.ima = None
+  def __init__(self,fname,alpha=1.,Nppr=64,ImgRescale=1,was_alpha_used_in_measurment=False,is_density_weighted=True,filter_kspace=False):
+    """
+      initiate paramters for reconstruction of concept
+      fname : path to twix file
+      optional:
+        alpha=1, scale the FOV 
+                 set larger than 1 will reconstruct a larger FOV, the Image matrix  
+                 remains the same (i.e. lower resolution)
+        was_alpha_used_in_measurment=False, set to True if sequence .ini file 
+                 applied alpha to the chosen FOV (ensure correct FOV size)
+        Nppr=64, default is 64, number RO points pere ring
+        ImgRescale=1 rescale the image matrix (i.e. resolution)
+                 set to 2 to double reconstructed matrix size, maintaining the FOV
+        is_density_weighted (True), this is unused, and instead auto detected
+        filter_kspace (False), reduces ringing, useful for unwighted scans
+    """
     self.kdata = None
     self.img = None
     self.hdr = None
@@ -46,8 +48,9 @@ class Recon:
     self.Noffset = 8 # number of datapoints to offset RO by
     self.was_alpha_used_in_measurment = was_alpha_used_in_measurment
     self.Nti = 2 # assume two temporal interleaves
-    self.is_dw = is_density_weighted # assume density weighting
+    self.is_dw = is_density_weighted # assume density weighting, not this is unused and instead auto detected
     self.filter_kspace = filter_kspace
+    self.sens = None
     
   """ 
     Read in raw data and header from twix file
@@ -201,7 +204,6 @@ class Recon:
         density_weighting = y*density_weighting
         #
     
-    
     kca = np.exp(1j*(np.linspace(0,ppr-1,ppr).reshape(ppr,1)*2*np.pi/ppr))* r.reshape(1,self.NRings)
     kca[:,::2] = kca[:,::2]*np.exp(1j*np.pi/ppr) # % shift odd rings by half a k point
     
@@ -213,7 +215,7 @@ class Recon:
   """
   Generate an image with these dimentions
   [x][y][z][TI,TIME][Cha]
-  todo swap time and cha
+  
   """
   def calcImg(self):
     if not isinstance( self.kdata,np.ndarray):
@@ -260,11 +262,68 @@ class Recon:
     # permute to make dimentions confirm to nifti with first three as spatial 
     # x,y,z,(3)time,(4)coil
     self.img = self.img.transpose([2,3,4,0,1])
-    #return self.img
+    self.ImgSz = self.img.shape
     
     return True
     
-  def saveIma2Nii(self,fname='',print_headers=True):
+  def loadSensMaps(self,fname_nii):
+    """
+    input:
+        fname_nii, file name to load sensitivity maps from
+    effect: 
+        sens maps loaded to self.sens
+    """
+    sens = np.asanyarray(nib.load(fname_nii).dataobj)
+    # test data consistent, X, Y and CH
+    if sens.shape == (self.ImgSz[0:3]+(1,self.ImgSz[4])):
+        self.sens = sens
+    else:
+        print('Cannot use sensitivies, dimentions do not agree, disgarding them')
+        print(sens.shape)
+ 
+    
+  def combineCoils(self):
+    """
+        if no maps already loaded, calculate maps using first data point
+        it is strongly recomended to use sensitvity maps calculated from a water spectra
+    effect:
+        sens maps recorded in self.sens
+        images combined with self.sens
+    """
+    if isinstance(self.sens, type(None)):
+        print('Calculating coil sensitivies with espirit')
+        sens = espirit.espirit_img(self.img[:,:,:,0,:].transpose((3,0,1,2)),dims=self.ImgSz[0:2])
+        self.sens = sens.reshape([self.ImgSz[4]]+list(self.ImgSz[0:3])+[1]).transpose((1,2,3,4,0))
+  
+    self.img = np.sum(self.img*np.conj(self.sens),axis=4)
+    
+  def saveSensMaps(self,fname_nii=''):
+    """
+    input:
+        fname_nii, file name to save sensitivity maps to, if none given the raw data file appended with '_sens' will be used
+    
+    effect: 
+        file saved
+    """
+    if len(fname_nii)==0:
+        fname_nii = self.fname_raw.replace('.dat','_sens.nii.gz')
+        
+    if isinstance(self.sens, type(None)):
+       print('No sens maps generated, skipping save')
+    else: 
+      self.saveIma2Nii(fname=fname_nii,print_headers=False,img_to_save=self.sens)
+  
+  def saveIma2Nii(self,fname='',print_headers=True,img_to_save=None):
+    """
+    input:
+        fname, file name to save to nifti, if none given the rawdata filename with .nii.gz will be used
+        print_headers: True to print info saved to nii
+        img_to_save: used to save alternative images to an nii file such as sens maps
+    effect: 
+        file saved
+    """
+    if isinstance(img_to_save, type(None)):
+        img_to_save = self.img.conj()
     
     p = self.hdr.Phoenix
     m = self.hdr.Meas
@@ -346,7 +405,7 @@ class Recon:
     
     # Calculate the Orientation
     qform = self.getOrientation()
-    newobj = nib.nifti2.Nifti2Image(self.img.conj(), qform)
+    newobj = nib.nifti2.Nifti2Image(img_to_save, qform)
     
     # Set q_form >0
     newobj.header.set_qform(qform)
@@ -368,15 +427,19 @@ class Recon:
     # # From nii obj and write    
     nib.save(newobj,fname)
 
-  """
-  Return qform matrix for nifit (4x4 affine transformation form pixel to patient space)
-  The convention is: +x = Right; +y = Anterior; +z = Superior. That is, moving in the positive x-direction (so that the x-coordinate increases) will be moving to the right, et
-  Siemens coordinates are SAG (Left to Right), COR (Anterior to Posterior), TRA (Foot to Head)
-  """
+
   def getOrientation(self):
+    """
+    output:
+      qform: matrix for nifit header (4x4 affine transformation form pixel to patient space)
+    
+    notes:
+    The convention is: +x = Right; +y = Anterior; +z = Superior. That is, moving in the positive x-direction (so that the x-coordinate increases) will be moving to the right, et
+    Siemens coordinates are SAG (Left to Right), COR (Anterior to Posterior), TRA (Foot to Head)
+    """
     phs,read,norm,pos = siemens.getOrientationVectors(self.hdr)
     qform = np.zeros([4,4])
-    res = self.FOV/self.ImgSz[-2]
+    res = self.FOV/self.ImgSz[0]
     # must point to center of voxel (i.e. need 
     fov_shift_phs = np.array(phs)*(self.FOV-res)/2
     fov_shift_read =  np.array(read)*(self.FOV-res)/2
@@ -393,6 +456,10 @@ class Recon:
     return qform
 
   def saveIma2mat(self,fname=''):
+    """
+    effect:
+    Save self.img to matlab variable (.mat) using file name or raw data file.mat
+    """
     from scipy import io
     if len(fname)==0:
       fname = self.fname_raw.replace('.dat','.mat')
@@ -406,6 +473,22 @@ class Recon:
 """
 
 def recon(fname_in,fname_out='',**kwargs):
+  """
+  Assume all default parameters, perform all steps
+  input: 
+    fname_in: name of raw data file to use
+  input (optional):
+    fname_out:  name of nifti file to save
+    alpha=1, scale the FOV 
+             set larger than 1 will reconstruct a larger FOV, the Image matrix  
+             remains the same (i.e. lower resolution)
+    was_alpha_used_in_measurment=False, set to True if sequence .ini file 
+             applied alpha to the chosen FOV (ensure correct FOV size)
+    Nppr=64, default is 64, number RO points pere ring
+    ImgRescale=1 rescale the image matrix (i.e. resolution)
+             set to 2 to double reconstructed matrix size, maintaining the FOV
+    filter_kspace (False), this reduces, useful for unwighted scans
+  """
   my_recon = Recon(fname_in,**kwargs)
   my_recon.getKdata()
   my_recon.arrangeData()
@@ -414,7 +497,7 @@ def recon(fname_in,fname_out='',**kwargs):
   return my_recon
   
 def recon_uw(fname_in,fname_out='',**kwargs):
-  my_recon = Recon(fname_in,is_density_weighted=False,**kwargs)
+  my_recon = Recon(fname_in,is_density_weighted=False,filter_kspace=True,**kwargs)
   my_recon.getKdata()
   my_recon.arrangeData()
   my_recon.calcImg()

espirit.py

+# Mark Chiew
+# mark.chiew@ndcn.ox.ac.uk
+import numpy as np
+
+fft  = lambda x, ax : np.fft.fftshift(np.fft.fftn(np.fft.ifftshift(x, axes=ax), axes=ax, norm='ortho'), axes=ax) 
+ifft = lambda X, ax : np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(X, axes=ax), axes=ax, norm='ortho'), axes=ax) 
+
+
+def espirit(x, dims=(128,128), kernel=(5,5), eig_thresh=0.02, mask_thresh=0.99):
+    """
+    Input:
+        x: [Nc, Kx, Ky, Nz]
+        dims: (Nx, Ny)
+
+    Output:
+        sens: [Nx, Ny, Nz, Nc]
+    """
+
+    if x.ndim < 4:
+        x = x[:,:,:,np.newaxis]
+    
+    # Get dimensions
+    Nc = x.shape[0]
+    Nz = x.shape[3]
+    Kx = x.shape[1]-kernel[0]+1
+    Ky = x.shape[2]-kernel[1]+1
+    pad_x = [np.int(i) for i in (np.ceil((dims[0]-kernel[0])/2), np.floor((dims[0]-kernel[0])/2))]
+    pad_y = [np.int(i) for i in (np.ceil((dims[1]-kernel[1])/2), np.floor((dims[1]-kernel[1])/2))]
+    sens  = np.zeros((Nc, dims[0], dims[1], Nz), dtype='complex')
+    m  = np.zeros((dims[0], dims[1], Nz), dtype='complex')
+
+    # Loop over slices (Nz)
+    for z in range(Nz):
+        # Initialise Hankel matrix
+        H = np.zeros((Nc, np.prod(kernel), Kx*Ky), dtype='complex')
+        for i in range(Kx):
+            for j in range(Ky):
+                # Populate Hankel matrix
+                H[:,:,i*Ky+j] = x[:,i:i+kernel[0],j:j+kernel[1],z].reshape((Nc,-1))
+
+        # Only keep svd singular vectors corresponding to above threshold singular values
+        U,S,_ = np.linalg.svd(H.reshape((-1,Kx*Ky)), full_matrices=False)
+        U = U[:,S>S[0]*eig_thresh]   
+
+        # Get zero-paded kernel images
+        U = U.reshape((Nc,kernel[0],kernel[1],-1))
+        U = np.pad(U, ((0,0), pad_x, pad_y, (0,0)))
+        U = np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(U, axes=(1,2)), axes=(1,2)), axes=(1,2))
+        U = U*np.prod(dims)/np.sqrt(np.prod(kernel))
+
+        # Perform voxel-wise SVD on coil x component matrices, voxelwise, keeping 1st component
+        for i in range(dims[0]):
+            for j in range(dims[1]):
+                W,S,_ = np.linalg.svd(U[:,i,j,:], full_matrices=False)
+                sens[:,i,j,z] = W[:,0]
+                m[i,j,z] = S[0]
+     # rotate phase relative to first channel and return
+    return np.squeeze(sens*np.exp(-1j*np.angle(sens[[0],:,:,:]))*(np.abs(m)>mask_thresh))
+    
+    
+ 
+def espirit_img(ximg, dims=(128,128), kernel=(5,5), eig_thresh=0.02, mask_thresh=0.99):
+    """
+    return sensity maps with images as input
+    Input:
+        ximg: [Nc, x, y, z]
+        dims: (Nx, Ny)
+
+    Output:
+        sens: [Nx, Ny, Nz, Nc]
+    """
+    return espirit(fft(ximg,(1,2,3)),dims,kernel,eig_thresh,mask_thresh)
\ No newline at end of file