Source code for mri.operators.fourier.utils.processing
# -*- coding: utf-8 -*-
##########################################################################
# pySAP - Copyright (C) CEA, 2017 - 2018
# Distributed under the terms of the CeCILL-B license, as published by
# the CEA-CNRS-INRIA. Refer to the LICENSE file or to
# http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
# for details.
##########################################################################
"""
Common tools for MRI image reconstruction.
"""
# System import
import warnings
# Third party import
import numpy as np
import scipy.fftpack as pfft
from scipy.interpolate import griddata, RegularGridInterpolator
from mrinufft import get_density
[docs]def convert_mask_to_locations(mask):
"""Return the converted Cartesian mask as sampling locations.
Parameters
----------
mask: numpy.ndarray, {0,1}
ND matrix, not necessarly a square matrix.
Returns
-------
samples_locations: numpy.ndarray
samples location between [-0.5, 0.5[ of shape MxN where M is the
number of 1 values in the mask.
"""
locations = np.where(mask == 1)
rslt = []
for dim, loc in enumerate(locations):
loc_n = loc.astype("float") / mask.shape[dim] - 0.5
rslt.append(loc_n)
return np.asarray(rslt).T
[docs]def convert_locations_to_mask(samples_locations, img_shape):
"""Return the converted the sampling locations as Cartesian mask.
Parameters
----------
samples_locations: numpy.ndarray
samples locations between [-0.5, 0.5[.
img_shape: tuple of int
shape of the desired mask, not necessarly a square matrix.
Returns
-------
mask: numpy.ndarray, {0,1}
2D matrix, not necessarly a square matrix.
"""
if samples_locations.shape[-1] != len(img_shape):
raise ValueError("Samples locations dimension doesn't correspond to ",
"the dimension of the image shape")
locations = np.copy(samples_locations).astype("float")
locations += 0.5
locations = np.clip(locations, 0, 0.999)
locations *= img_shape
locations = locations.astype("int")
mask = np.zeros(img_shape, dtype="int")
mask[tuple(locations.T)] = 1
return mask
[docs]def normalize_frequency_locations(samples, Kmax=None):
"""
This function normalizes the sample locations between [-0.5; 0.5[ for
the non-Cartesian case.
Parameters
----------
samples: numpy.ndarray
Unnormalized samples
Kmax: int, float, array-like or None
Maximum Frequency of the samples locations is supposed to be equal to
base Resolution / (2* Field of View)
Returns
-------
normalized_samples: numpy.ndarray
Same shape as the parameters but with values between [-0.5; 0.5[
"""
samples_locations = np.copy(samples.astype('float'))
if Kmax is None:
Kmax = np.abs(samples_locations).max(axis=0)
elif isinstance(Kmax, (float, int)):
Kmax = [Kmax] * samples_locations.shape[-1]
Kmax = np.array(Kmax)
samples_locations /= (2 * Kmax)
if np.abs(samples_locations).max() >= 0.5:
warnings.warn("Frequencies outside the 0.5 limit.")
return samples_locations
[docs]def discard_frequency_outliers(kspace_loc, kspace_data):
"""
This function discards the samples outside [-0.5; 0.5[ for
the non-Cartesian case.
Parameters
----------
kspace_loc: numpy.ndarray
The sample locations previously normalized around [-0.5; 0.5[
using Kmax.
kspace_data: numpy.ndarray
The samples corresponding to kspace_loc defined above.
Returns
-------
reduced_kspace_loc: numpy.ndarray
The sample locations reduced strictly to [-0.5; 0.5[ by discarding
outliers.
reduced_kspace_data: numpy.ndarray
The samples corresponding to reduced_kspace_loc defined above.
"""
kspace_mask = np.all((kspace_loc < 0.5) & (kspace_loc >= -0.5), axis=-1)
kspace_loc = kspace_loc[kspace_mask]
kspace_data = kspace_data[:, kspace_mask]
return np.ascontiguousarray(kspace_loc), np.ascontiguousarray(kspace_data)
[docs]def gridded_inverse_fourier_transform_nd(kspace_loc,
kspace_data, grid, method):
"""
This function calculates the gridded Inverse fourier transform
from Interpolated non-Cartesian data into a cartesian grid
Parameters
----------
kspace_loc: numpy.ndarray
The N-D k_space locations of size [M, N]
kspace_data: numpy.ndarray
The k-space data corresponding to k-space_loc above
grid: numpy.ndarray
The Gridded matrix for which you want to calculate k_space Smaps
method: {'linear', 'nearest', 'cubic'}
Method of interpolation for more details see scipy.interpolate.griddata
documentation
Returns
-------
np.ndarray
The gridded inverse fourier transform of given kspace data
"""
gridded_kspace = griddata(kspace_loc,
kspace_data,
grid,
method=method,
fill_value=0)
return np.swapaxes(pfft.fftshift(
pfft.ifftn(pfft.ifftshift(gridded_kspace))), 1, 0)
[docs]def gridded_inverse_fourier_transform_stack(kspace_data_sorted,
kspace_plane_loc,
idx_mask_z,
grid,
volume_shape,
method):
"""
This function calculates the gridded Inverse fourier transform
from Interpolated non-Cartesian data into a cartesian grid. However,
the IFFT is done similar to Stacked Fourier transform.
We expect the kspace data to be limited to a grid on z, we calculate
the inverse fourier transform by-
1) Grid data in each plane (for all points in a plane)
2) Interpolate data along z, if we have undersampled data along z
3) Apply an IFFT on the 3D data that was gridded and interpolated in z.
Parameters
----------
kspace_data_sorted: numpy.ndarray
The sorted k-space data corresponding to kspace_plane_loc above
kspace_plane_loc: numpy.ndarray
The N-D k_space locations of size [M, N]. These hold locations only
in plane, extracted using get_stacks_fourier function
idx_mask_z: numpy.ndarray
contains the indices of the acquired Fourier plane. Extracted using
get_stacks_fourier function
grid: tuple
The Gridded matrix for which you want to calculate k_space Smaps.
Should be given as a tuple of ndarray
volume_shape: tuple
Reconstructed volume shape
method: {'linear', 'nearest', 'cubic'}, optional
Method of interpolation for more details see scipy.interpolate.griddata
documentation
Returns
-------
np.ndarray
The gridded inverse fourier transform of given kspace data
"""
gridded_kspace = np.zeros(volume_shape, dtype=kspace_data_sorted.dtype)
stack_len = len(kspace_plane_loc)
for i, idx_z in enumerate(idx_mask_z):
gridded_kspace[:, :, idx_z] = griddata(
kspace_plane_loc,
kspace_data_sorted[i*stack_len:(i+1)*stack_len],
grid,
method=method,
fill_value=0,
)
# Check if we have undersampled in Z direction, in which case,
# we need to interpolate values along z to get a good reconstruction.
if len(idx_mask_z) < volume_shape[2]:
# Interpolate along z direction
grid_loc = [
np.linspace(-0.5, 0.5, volume_shape[i], endpoint=False)
for i in range(3)
]
interp_z = RegularGridInterpolator(
(*grid_loc[0:2], grid_loc[2][idx_mask_z]),
gridded_kspace[:, :, idx_mask_z],
bounds_error=False,
fill_value=None,
)
unsampled_z = list(
set(np.arange(volume_shape[2])) - set(idx_mask_z)
)
mask = np.zeros(volume_shape)
mask[:, :, unsampled_z] = 1
loc = convert_mask_to_locations(mask)
gridded_kspace[:, :, unsampled_z] = np.reshape(
interp_z(loc),
(*volume_shape[0:2], len(unsampled_z)),
)
# Transpose every image in each slice
return np.swapaxes(np.fft.fftshift(np.fft.ifftn(np.fft.ifftshift(
gridded_kspace))), 0, 1)
[docs]def check_if_fourier_op_uses_sense(fourier_op):
"""Utils function to check if fourier operator uses SENSE recon
Parameters
----------
fourier_op: object of class FFT, NonCartesianFFT or Stacked3DNFFT in
mri.operators
the fourier operator for which we want to check if SENSE is
supported
Returns
-------
bool
True if SENSE recon is being used
"""
from ..non_cartesian import NonCartesianFFT
if isinstance(fourier_op, NonCartesianFFT):
return fourier_op.impl.uses_sense
else:
return False
[docs]def estimate_density_compensation(kspace_loc, volume_shape, implementation='pipe', **kwargs):
""" Utils function to obtain the density compensator for a
given set of kspace locations.
Parameters
----------
kspace_loc: numpy.ndarray
the kspace locations
volume_shape: numpy.ndarray
the volume shape
implementation: str default 'pipe'
the implementation of the non-cartesian operator
can be 'pipe' which needs gpuNUFFT or 'cell_count'
kwargs: dict
extra keyword arguments to be passed to the density
compensation estimation
"""
density_comp = get_density(implementation)(
kspace_loc,
volume_shape,
**kwargs
)
return np.squeeze(density_comp)