mfield.py

from __future__ import division

__all__ = ['genDBz', 'genDf', 'genSusMap', 'gamma2pi']

import numpy as np
from scipy.special import erf


gamma2pi = 42.57747892e6 # gyromagentic ratio
chiFerritin = 1.30e-9 # Ferritin susceptibility following Schenck 96 
                      # (*e3 because we hace [c] in \mu g / g and not mg / g )


# Function that returns the B-field change from a given susceptibility map
# Checked against Marques & Bowtell paper
def genDBz(susMap, B0, axis, phi = 0, theta = np.pi/2):
    """ 
    Generates magnetic field shift resulting from action of magnetic field B0 along axis specified by axis and/or angles phi and theta.
    susMap: 3D array of susceptibility in SI units
    B0: magnetic field in Tesla
    axis: orientation of B0 field along x(axis=0)-, y(axis=1)- or z(axis=2)-axis; or along arbitrary axis(=3) specified by polar angle theta and azimuthal angle phi w.r.t. z-axis.
    phi: if axis==3: azimuthal angle of B0
    theta: if axis==3: polar angle of B0
    """
    # Coordinate grid for kernel calculation
    x = np.linspace(-susMap.shape[0]+1, susMap.shape[0]-1, susMap.shape[0])
    y = np.linspace(-susMap.shape[1]+1, susMap.shape[1]-1, susMap.shape[1])
    z = np.linspace(-susMap.shape[2]+1, susMap.shape[2]-1, susMap.shape[2])
    
    kx, ky, kz = np.meshgrid(x,y,z)
    
    if axis == 0:
        kernel = 1/3. - kx*kx / (kx*kx + ky*ky + kz*kz) 
    elif axis == 1:
        kernel = 1/3. - ky*ky / (kx*kx + ky*ky + kz*kz) 
    elif axis == 2:
        kernel = 1/3. - kz*kz / (kx*kx + ky*ky + kz*kz) 
    elif axis == 3: # kernel for B field with polar angle theta and azimuthal angle phi
        kernel = 1/3. - ((kx*np.cos(phi) + ky*np.sin(phi))*np.sin(theta) + kz*np.cos(theta))**2 / (kx*kx + ky*ky + kz*kz) 
        
    kernel[susMap.shape[0]//2,susMap.shape[1]//2,susMap.shape[2]//2] = 0 # exclude sphere of Lorentz 

    mu0Mz = susMap * B0
    Fmu0Mz = np.fft.fftshift(np.fft.fftn(np.fft.fftshift(mu0Mz)))

    kFmu0Mz = Fmu0Mz * kernel # = B~ dz (k) is FT of mu0Mz with dipole kernel, [6] Bowtell
    kFmu0Mz[susMap.shape[0]//2,susMap.shape[1]//2,susMap.shape[2]//2] = B0*susMap.mean()\
                                                                        /3.*susMap.shape[0]*susMap.shape[1]*susMap.shape[2]
        # following Bowtell => sphere of Lorentz at the zero frequency position after ifftshift, reflects mean field strength
    return np.real(np.fft.ifftshift(np.fft.ifftn(np.fft.ifftshift(kFmu0Mz)))) #Bz by fourier backtrafo

def genDf(susMap, B0, axis, phi = 0, theta = np.pi/2):
    """ 
    Generates local Larmor frequency shift resulting from action of magnetic field B0 along axis specified by axis and/or angles phi and theta.
    susMap: 3D array of susceptibility in SI units
    B0: magnetic field in Tesla
    axis: orientation of B0 field along x(axis=0)-, y(axis=1)- or z(axis=2)-axis; or along arbitrary axis specified by polar angle theta and azimuthal angle phi w.r.t. z-axis.
    phi: if axis==3: azimuthal angle of B0
    theta: if axis==3: polar angle of B0
    """
    gamma2pi = 42.57747892e6 # gyromagentic ratio

    return genDBz(susMap, B0, axis,phi, theta)*gamma2pi # f = gamma / 2pi * dB

def phi(x): #needed for smoothing
    return 1/2.*(1 + erf(x/(2*np.sqrt(2))))


def genSusMap(cMap, sliceHeight,stackedaxis, drytissuecorrection = False, smooth = False):
    # Creates a correctly spaced volume of the samples susceptibility assuming that it is proportional to the 
    # iron density (like Ferritin).
    drytissueprop = 0.2 # Proportion of dry tissue in physiological conditions
    
    if drytissuecorrection:
        chiFerritin = chiFerritin*drytissueprop
    
    suscep = cMap * mfield.chiFerritin

    # To real scaled cube: sliceHeight high slices with area of cMap.shape[1]*cMap.shape[2] px
    if stackedaxis==0:
        suscep3dz = np.zeros([sliceHeight * cMap.shape[0],cMap.shape[1],cMap.shape[2]])
        for i in range(cMap.shape[0]):
            suscep3dz[(i*sliceHeight):((i+1)*sliceHeight),:,:] = suscep[i,:,:]

    elif stackedaxis==1:
        suscep3dz = np.zeros([cMap.shape[0],sliceHeight * cMap.shape[1],cMap.shape[2]])
        for i in range(cMap.shape[1]):
            suscep3dz[:,(i*sliceHeight):((i+1)*sliceHeight),:] = suscep[:,i,:]

    elif stackedaxis==2:
        suscep3dz = np.zeros([cMap.shape[0],cMap.shape[1],sliceHeight * cMap.shape[2]])
        for i in range(cMap.shape[2]):
            suscep3dz[:,:,(i*sliceHeight):((i+1)*sliceHeight)] = suscep[:,:,i][:,:,None]
    
    if smooth and sliceHeight==4 and cMap.shape[0]==10 and stackedaxis==2: # For "original" 10 mu slices.
        for i in range(10):
            susMap[:,:,4*i-1]= phi(1)*susMap[:,:,4*i-1] + phi(-1)*susMap[:,:,4*i]
            susMap[:,:,4*i] =  phi(-1)*susMap[:,:,4*i-1] + phi(1)*susMap[:,:,4*i]
        print('sucessfully smoothed')
    else:
        print('smoothing was not done')
    return suscep3dz