Source code for mrinufft.io.cfl

"""utility functions for Complex-Float data file.

This allow compatibility with BART toolbox and others.

References
----------
BART Toolbox: https://bart-doc.readthedocs.io/en/latest/data.html
"""

from pathlib import Path
import os
import mmap
import warnings
import numpy as np


def _readcfl(cfl_file, hdr_file=None):
    """Read a pair of .cfl/.hdr file to get a complex numpy array.

    Adapted from the BART python cfl library.

    Parameters
    ----------
    name : str
        Name of the file to read.

    Returns
    -------
    array : array_like
        Complex array read from the file.
    """
    basename = Path(cfl_file).with_suffix("")
    if hdr_file is None:
        hdr_file = basename.with_suffix(".hdr")
    cfl_file = basename.with_suffix(".cfl")

    # get dims from .hdr
    with open(hdr_file) as h:
        h.readline()  # skip
        line = h.readline()
    dims = [int(i) for i in line.split()]

    # remove singleton dimensions from the end
    n = np.prod(dims)
    dims_prod = np.cumprod(dims)
    dims = dims[: np.searchsorted(dims_prod, n) + 1]

    # load data and reshape into dims
    with open(cfl_file, "rb") as d:
        a = np.fromfile(d, dtype=np.complex64, count=n)
    return a.reshape(dims, order="F")


def _writecfl(array, cfl_file, hdr_file=None):
    """Write a pair of .cfl/.hdr file representing a complex array.

    Adapted from the BART python cfl library.

    Parameters
    ----------
    name : str
        Name of the file to write.
    array : array_like
        Array to write to file.

    """
    basename = Path(cfl_file).with_suffix("")
    if hdr_file is None:
        hdr_file = basename.with_suffix(".hdr")
    cfl_file = basename.with_suffix(".cfl")

    with open(hdr_file, "w") as h:
        h.write("# Dimensions\n")
        for i in array.shape:
            h.write("%d " % i)
        h.write("\n")

    size = np.prod(array.shape) * np.dtype(np.complex64).itemsize

    with open(cfl_file, "w+b") as d:
        os.ftruncate(d.fileno(), size)
        mm = mmap.mmap(d.fileno(), size, flags=mmap.MAP_SHARED, prot=mmap.PROT_WRITE)
        array = array.astype(np.complex64, copy=False)
        mm.write(np.ascontiguousarray(array.T))
        mm.close()


[docs] def traj2cfl(traj, shape, basename): """ Export a trajectory defined in MRI-nufft to a BART compatible format. The trajectory will be normalized to -(FOV-1)/2 +(FOV-1)/2, and reshape to BART format. Parameters ---------- traj: array_like trajectory array, shape (N_shot, N_points, 2 or 3) shape: tuple volume shape (FOV) """ traj_ = traj * (np.array(shape) - 1) if traj.shape[-1] == 2: traj_3d = np.zeros(traj_.shape[:-1] + (3,), dtype=traj_.dtype) traj_3d[..., :2] = traj_ traj_ = traj_3d else: traj_ = traj_.astype(np.complex64, copy=False) traj_ = traj_[None, None, ...] _writecfl(traj_.T, basename)
[docs] def cfl2traj(basename, shape=None): """Convert a trajectory BART file to a numpy array compatible with MRI-nufft. Parameters ---------- filename: str Base filename for the trajectory shape: optional Shape of the Image domain. Returns ------- np.ndarray MRI-NUFFT compatible trajectory of shape (n_shot, n_samples, dim) """ traj_raw = _readcfl(basename) # Convert to float array and take only the real part traj = np.ascontiguousarray(traj_raw.T.view("(2,)float32")[..., 0]) if np.all(traj[..., -1] == 0): warnings.warn("2D Trajectory Detected") traj = traj[..., :-1] if shape is None: maxs = [np.max(traj[..., i]) for i in range(traj.shape[-1])] mins = [np.min(traj[..., i]) for i in range(traj.shape[-1])] shape = np.array(maxs) - np.array(mins) warnings.warn(f"Estimated shape {shape}") else: shape = np.asarray(shape) - 1 traj /= np.asarray(shape) return np.squeeze(traj)