"""
NUFFT CPU class (deprecated)
=======================================
"""
from __future__ import absolute_import
import numpy
import scipy.sparse
import numpy.fft
import scipy.linalg
import scipy.special
from ..src._helper import helper, helper1
class NUFFT_cpu:
"""
Class NUFFT_cpu
"""
def __init__(self):
"""
Constructor.
:param None:
:type None: Python NoneType
:return: NUFFT: the pynufft_hsa.NUFFT instance
:rtype: NUFFT: the pynufft_hsa.NUFFT class
:Example:
>>> from pynufft import NUFFT_cpu
>>> NufftObj = NUFFT_cpu()
"""
self.dtype = numpy.complex64 # : initial value: numpy.complex64
self.debug = 0 #: initial value: 0
self.Nd = () # : initial value: ()
self.Kd = () # : initial value: ()
self.Jd = () #: initial value: ()
self.ndims = 0 # : initial value: 0
self.ft_axes = () # : initial value: ()
self.batch = None # : initial value: None
pass
def plan(self, om, Nd, Kd, Jd, ft_axes=None):
"""
Plan the NUFFT_cpu object with the geometry provided.
:param om: The M off-grid locates in the frequency domain,
which is normalized between [-pi, pi]
:param Nd: The matrix size of the equispaced image.
Example: Nd=(256,256) for a 2D image;
Nd = (128,128,128) for a 3D image
:param Kd: The matrix size of the oversampled frequency grid.
Example: Kd=(512,512) for 2D image;
Kd = (256,256,256) for a 3D image
:param Jd: The interpolator size.
Example: Jd=(6,6) for 2D image;
Jd = (6,6,6) for a 3D image
:param ft_axes: (Optional) The axes for Fourier transform.
The default is all axes if 'None' is given.
:param batch: (Optional) Batch mode.
If the batch is provided, the last appended axis is the number
of identical NUFFT to be transformed.
The default is 'None'.
:type om: numpy.float array, matrix size = M * ndims
:type Nd: tuple, ndims integer elements.
:type Kd: tuple, ndims integer elements.
:type Jd: tuple, ndims integer elements.
:type ft_axes: None, or tuple with optional integer elements.
:type batch: None, or integer
:returns: 0
:rtype: int, float
:ivar Nd: initial value: Nd
:ivar Kd: initial value: Kd
:ivar Jd: initial value: Jd
:ivar ft_axes: initial value: None
:ivar batch: initial value: None
:Example:
>>> from pynufft import NUFFT_cpu
>>> NufftObj = NUFFT_cpu()
>>> NufftObj.plan(om, Nd, Kd, Jd)
or
>>> NufftObj.plan(om, Nd, Kd, Jd, ft_axes, batch)
"""
self.ndims = len(Nd) # : initial value: len(Nd)
if ft_axes is None:
ft_axes = range(0, self.ndims)
self.ft_axes = ft_axes # default: all axes (range(0, self.ndims)
self.st = helper.plan(om, Nd, Kd, Jd, ft_axes=ft_axes,
format='CSR')
self.Nd = self.st['Nd'] # backup
self.Kd = self.st['Kd']
# backup
self.sn = numpy.asarray(self.st['sn'].astype(self.dtype), order='C')
if batch is None: # single-coil
self.parallel_flag = 0
self.batch = 1
else: # multi-coil
self.parallel_flag = 1
self.batch = batch
if self.parallel_flag == 1:
self.multi_Nd = self.Nd + (self.batch, )
self.uni_Nd = self.Nd + (1, )
self.multi_Kd = self.Kd + (self.batch, )
self.multi_M = (self.st['M'], ) + (self.batch, )
self.multi_prodKd = (numpy.prod(self.Kd), self.batch)
self.sn = numpy.reshape(self.sn, self.Nd + (1,))
elif self.parallel_flag == 0:
self.multi_Nd = self.Nd # + (self.Reps, )
self.uni_Nd = self.Nd
self.multi_Kd = self.Kd # + (self.Reps, )
self.multi_M = (self.st['M'], )
self.multi_prodKd = (numpy.prod(self.Kd), )
# Calculate the density compensation function
self.sp = self.st['p'].copy().tocsr()
self.spH = (self.st['p'].getH().copy()).tocsr()
self.Kdprod = numpy.int32(numpy.prod(self.st['Kd']))
self.Jdprod = numpy.int32(numpy.prod(self.st['Jd']))
del self.st['p'], self.st['sn']
self.NdCPUorder, self.KdCPUorder, self.nelem = helper.preindex_copy(
self.st['Nd'],
self.st['Kd'])
self.volume = {}
self.volume['cpu_coil_profile'] = numpy.ones(self.multi_Nd)
return 0
def _precompute_sp(self):
"""
Private: Precompute adjoint (gridding) and Toepitz interpolation
matrix.
:param None:
:type None: Python Nonetype
:return: self: instance
"""
try:
W0 = numpy.ones((self.st['M'],), dtype=numpy.complex64)
W = self.xx2k(self.adjoint(W0))
self.W = (W*W.conj())**0.5
del W0
del W
except:
print("errors occur in self.precompute_sp()")
raise
def reset_sense(self):
self.volume['cpu_coil_profile'].fill(1.0)
def set_sense(self, coil_profile):
self.volume = {}
if coil_profile.shape == self.Nd + (self.batch, ):
self.volume['cpu_coil_profile'] = coil_profile
else:
print('The shape of coil_profile might be wrong')
print('coil_profile.shape = ', coil_profile.shape)
print('shape of Nd + (batch, ) = ', self.Nd + (self.batch, ))
def forward_one2many(self, x):
"""
Assume x.shape = self.Nd
"""
# try:
x2 = x.reshape(self.uni_Nd, order='C')*self.volume['cpu_coil_profile']
# except:
# x2 = x
y2 = self.forward(x2)
return y2
def adjoint_many2one(self, y):
"""
Assume y.shape = self.multi_M
"""
x2 = self.adjoint(y)
x = x2*self.volume['cpu_coil_profile'].conj()
try:
x3 = numpy.mean(x, axis=self.ndims)
except:
x3 = x
del x
return x3
def solve(self, y, solver=None, *args, **kwargs):
"""
Solve NUFFT_cpu.
:param y: data, numpy.complex64. The shape = (M,) or (M, batch)
:param solver: 'cg', 'L1TVOLS', 'lsmr', 'lsqr', 'dc', 'bicg',
'bicgstab', 'cg', 'gmres','lgmres'
:param maxiter: the number of iterations
:type y: numpy array, dtype = numpy.complex64
:type solver: string
:type maxiter: int
:return: numpy array with size.
The shape = Nd ('L1TVOLS') or Nd + (batch,)
('lsmr', 'lsqr', 'dc','bicg','bicgstab','cg', 'gmres','lgmres')
"""
from ..linalg.solve_cpu import solve
x2 = solve(self, y, solver, *args, **kwargs)
return x2 # solve(self, y, solver, *args, **kwargs)
def forward(self, x):
"""
Forward NUFFT on CPU
:param x: The input numpy array, with the size of Nd or Nd + (batch,)
:type: numpy array with the dtype of numpy.complex64
:return: y: The output numpy array, with the size of (M,) or (M, batch)
:rtype: numpy array with the dtype of numpy.complex64
"""
y = self.k2y(self.xx2k(self.x2xx(x)))
return y
def adjoint(self, y):
"""
Adjoint NUFFT on CPU
:param y: The input numpy array, with the size of (M,) or (M, batch)
:type: numpy array with the dtype of numpy.complex64
:return: x: The output numpy array,
with the size of Nd or Nd + (batch, )
:rtype: numpy array with the dtype of numpy.complex64
"""
x = self.xx2x(self.k2xx(self.y2k(y)))
return x
def selfadjoint_one2many2one(self, x):
y2 = self.forward_one2many(x)
x2 = self.adjoint_many2one(y2)
del y2
return x2
def selfadjoint(self, x):
"""
selfadjoint NUFFT on CPU
:param x: The input numpy array, with size=Nd
:type: numpy array with dtype =numpy.complex64
:return: x: The output numpy array, with size=Nd
:rtype: numpy array with dtype =numpy.complex64
"""
# x2 = self.adjoint(self.forward(x))
x2 = self.xx2x(self.k2xx(self.k2y2k(self.xx2k(self.x2xx(x)))))
return x2
def selfadjoint2(self, x):
try:
x2 = self.k2xx(self.W * self.xx2k(x))
except:
self._precompute_sp()
x2 = self.k2xx(self.W * self.xx2k(x))
return x2
def x2xx(self, x):
"""
Private: Scaling on CPU
Inplace multiplication of self.x_Nd by the scaling factor self.sn.
"""
xx = x * self.sn
return xx
def xx2k(self, xx):
"""
Private: oversampled FFT on CPU
Firstly, zeroing the self.k_Kd array
Second, copy self.x_Nd array to self.k_Kd array by cSelect
Third, inplace FFT
"""
output_x = numpy.zeros(self.multi_Kd, dtype=self.dtype, order='C')
for bat in range(0, self.batch):
output_x.ravel()[self.KdCPUorder * self.batch + bat] = xx.ravel()[
self.NdCPUorder * self.batch + bat]
k = numpy.fft.fftn(output_x, axes=self.ft_axes)
return k
def xx2k_one2one(self, xx):
"""
Private: oversampled FFT on CPU
First, zeroing the self.k_Kd array
Second, copy self.x_Nd array to self.k_Kd array by cSelect
Third, inplace FFT
"""
output_x = numpy.zeros(self.st['Kd'], dtype=self.dtype, order='C')
# for bat in range(0, self.batch):
output_x.ravel()[self.KdCPUorder] = xx.ravel()[self.NdCPUorder]
k = numpy.fft.fftn(output_x, axes=self.ft_axes)
return k
def k2vec(self, k):
k_vec = numpy.reshape(k, self.multi_prodKd, order='C')
return k_vec
def vec2y(self, k_vec):
'''
gridding:
'''
y = self.sp.dot(k_vec)
# y = self.st['ell'].spmv(k_vec)
return y
def k2y(self, k):
"""
Private: interpolation by the Sparse Matrix-Vector Multiplication
"""
y = self.vec2y(self.k2vec(k))
# numpy.reshape(self.st['p'].dot(Xk), (self.st['M'], ), order='F')
return y
def y2vec(self, y):
'''
regridding non-uniform data (unsorted vector)
'''
# k_vec = self.st['p'].getH().dot(y)
k_vec = self.spH.dot(y)
# k_vec = self.st['ell'].spmvH(y)
return k_vec
def vec2k(self, k_vec):
'''
Sorting the vector to k-spectrum Kd array
'''
k = numpy.reshape(k_vec, self.multi_Kd, order='C')
return k
def y2k(self, y):
"""
Private: gridding by the Sparse Matrix-Vector Multiplication
"""
k_vec = self.y2vec(y)
k = self.vec2k(k_vec)
return k
def k2xx(self, k):
"""
Private: the inverse FFT and image cropping (which is the reverse of
_xx2k() method)
"""
# dd = numpy.size(self.Kd)
k = numpy.fft.ifftn(k, axes=self.ft_axes)
xx = numpy.zeros(self.multi_Nd, dtype=self.dtype, order='C')
for bat in range(0, self.batch):
xx.ravel()[self.NdCPUorder * self.batch + bat] = k.ravel()[
self.KdCPUorder * self.batch + bat]
# xx = xx[crop_slice_ind(self.Nd)]
return xx
def k2xx_one2one(self, k):
"""
Private: the inverse FFT and image cropping
(which is the reverse of _xx2k() method)
"""
# dd = numpy.size(self.Kd)
k = numpy.fft.ifftn(k, axes=self.ft_axes)
xx = numpy.zeros(self.st['Nd'], dtype=self.dtype, order='C')
# for bat in range(0, self.batch):
xx.ravel()[self.NdCPUorder] = k.ravel()[self.KdCPUorder]
# xx = xx[crop_slice_ind(self.Nd)]
return xx
def xx2x(self, xx):
"""
Private: rescaling, which is identical to the _x2xx() method
"""
x = self.x2xx(xx)
return x
def k2y2k(self, k):
"""
Private: the integrated interpolation-gridding by the Sparse
Matrix-Vector Multiplication
"""
Xk = self.k2vec(k)
# k = self.spHsp.dot(Xk)
# k = self.spH.dot(self.sp.dot(Xk))
k = self.y2vec(self.vec2y(Xk))
k = self.vec2k(k)
return k
