<commit_before>from __future__ import division, print_function, absolute_import

import numpy as np


def _image_tv(x, axis=0, n_points=3):
    """ Computes total variation (TV) of matrix x accross a given axis and
    along two directions.

    Parameters
    ----------
    x : 2D ndarray
        matrix x
    axis : int (0 or 1)
        Axis which TV will be calculated. Default a is set to 0.
    n_points : int
        Number of points to be included in TV calculation.

    Returns
    -------
    ptv : 2D ndarray
        Total variation calculated from the right neighbours of each point
    ntv : 2D ndarray
        Total variation calculated from the left neighbours of each point
    """
    xs = x.copy() if axis else x.T.copy()

    # Add copies of the data so that data extreme points are also analysed
    xs = np.concatenate((xs[:, (-n_points-1):], xs, xs[:, 0:(n_points+1)]),
                        axis=1)

    ptv = np.absolute(xs[:, (n_points+1):(-n_points-1)] -
                      xs[:, (n_points+2):(-n_points)])
    ntv = np.absolute(xs[:, (n_points+1):(-n_points-1)] -
                      xs[:, (n_points):(-n_points-2)])
    for n in range(1, n_points):
        ptv = ptv + np.absolute(xs[:, (n_points+1+n):(-n_points-1+n)] -
                                xs[:, (n_points+2+n):(-n_points+n)])
        ntv = ntv + np.absolute(xs[:, (n_points+1-n):(-n_points-1-n)] -
                                xs[:, (n_points-n):(-n_points-2-n)])

    if axis:
        return ptv, ntv
    else:
        return ptv.T, ntv.T


def _gibbs_removal_1d(x, axis=0, n_points=3):
    """ Suppresses Gibbs ringing along a given axis using fourier sub-shifts.

    Parameters
    ----------
    x : 2D ndarray
        Matrix x.
    axis : int (0 or 1)
        Axis in which Gibbs oscillations will be suppressed. Default is set
        to 0
    n_points : int, optional
        Number of neighbours to access local TV (see note). Default is set to
        3.

    Returns
    -------
    xc : 2D ndarray
        Matrix with suppressed Gibbs oscillations along the given axis.

    Note
    ----
    This function suppresses the effects of Gibbs oscillations based on the
    analysis of local total variation (TV). Although artefact correction is
    done based on two adjanced points for each voxel, total variation should be
    accessed in a larger range of neighbours. The number of the neighbours to
    be considered in TV calculation can be adjusted using parameter n_points.
    """
    ssamp = np.linspace(0.02, 0.9, num=45)

    xs = x.copy() if axis else x.T.copy()

    # TV for shift zero (baseline)
    tvr, tvl = _image_tv(xs, axis=1, n_points=n_points)
    tvp = np.minimum(tvr, tvl)
    tvn = tvp.copy()

    # Find optimal shift for gibbs removal
    isp = xs.copy()
    isn = xs.copy()
    sp = np.zeros(xs.shape)
    sn = np.zeros(xs.shape)
    N = xs.shape[1]
    c = np.fft.fftshift(np.fft.fft2(xs))
    k = np.linspace(-N/2, N/2-1, num=N)
    k = (2.0j * np.pi * k) / N
    for s in ssamp:
        # Access positive shift for given s
        img_p = abs(np.fft.ifft2(np.fft.fftshift(c * np.exp(k*s))))
        tvsr, tvsl = _image_tv(img_p, axis=1, n_points=n_points)
        tvs_p = np.minimum(tvsr, tvsl)

        # Access negative shift for given s
        img_n = abs(np.fft.ifft2(np.fft.fftshift(c * np.exp(-k*s))))
        tvsr, tvsl = _image_tv(img_n, axis=1, n_points=n_points)
        tvs_n = np.minimum(tvsr, tvsl)

        # Update positive shift params
        isp[tvp > tvs_p] = img_p[tvp > tvs_p]
        sp[tvp > tvs_p] = s
        tvp[tvp > tvs_p] = tvs_p[tvp > tvs_p]

        # Update negative shift params
        isn[tvn > tvs_n] = img_n[tvn > tvs_n]
        sn[tvn > tvs_n] = s
        tvn[tvn > tvs_n] = tvs_n[tvn > tvs_n]

    # check non-zero sub-voxel shifts
    idx = np.nonzero(sp + sn)

    # use positive and negative optimal sub-voxel shifts to interpolate to
    # original grid points
    xs[idx] = (isp[idx] - isn[idx])/(sp[idx] + sn[idx])*sn[idx] + isn[idx]

    return xs if axis else xs.T


def _weights(shape):
    """ Computes the weights necessary to combine two images processed by
    the 1D Gibbs removal procedure along two different axes [1]_.

    Parameters
    ----------
    shape : tuple
        shape of the image

    Returns
    -------
    G0 : 2D ndarray
        Weights for the image corrected along axis 0.
    G1 : 2D ndarray
        Weights for the image corrected along axis 1.

    References
    ----------
    .. [1] <NAME>, <NAME>, <NAME>, <NAME>. Gibbs-ringing artifact
           removal based on local subvoxel-shifts. Magn Reson Med. 2016
           doi: 10.1002/mrm.26054.
    """
    G0 = np.zeros(shape)
    G1 = np.zeros(shape)
    k0 = np.linspace(-np.pi, np.pi, num=shape[0])
    k1 = np.linspace(-np.pi, np.pi, num=shape[1])

    # Middle points
    K1, K0 = np.meshgrid(k1[1:-1], k0[1:-1])
    cosk0 = 1.0 + np.cos(K0)
    cosk1 = 1.0 + np.cos(K1)
    G1[1:-1, 1:-1] = cosk0 / (cosk0+cosk1)
    G0[1:-1, 1:-1] = cosk1 / (cosk0+cosk1)

    # Boundaries
    G1[1:-1, 0] = G1[1:-1, -1] = 1
    G1[0, 0] = G1[-1, -1] = G1[0, -1] = G1[-1, 0] = 1/2
    G0[0, 1:-1] = G0[-1, 1:-1] = 1
    G0[0, 0] = G0[-1, -1] = G0[0, -1] = G0[-1, 0] = 1/2

    return G0, G1


def _gibbs_removal_2d(image, n_points=3, G0=None, G1=None):<commit_msg>PEP8: minor pep8 fix in gibbs.py<commit_after>from __future__ import division, print_function, absolute_import

import numpy as np


def _image_tv(x, axis=0, n_points=3):
    """ Computes total variation (TV) of matrix x accross a given axis and
    along two directions.

    Parameters
    ----------
    x : 2D ndarray
        matrix x
    axis : int (0 or 1)
        Axis which TV will be calculated. Default a is set to 0.
    n_points : int
        Number of points to be included in TV calculation.

    Returns
    -------
    ptv : 2D ndarray
        Total variation calculated from the right neighbours of each point
    ntv : 2D ndarray
        Total variation calculated from the left neighbours of each point
    """
    xs = x.copy() if axis else x.T.copy()

    # Add copies of the data so that data extreme points are also analysed
    xs = np.concatenate((xs[:, (-n_points-1):], xs, xs[:, 0:(n_points+1)]),
                        axis=1)

    ptv = np.absolute(xs[:, (n_points+1):(-n_points-1)] -
                      xs[:, (n_points+2):(-n_points)])
    ntv = np.absolute(xs[:, (n_points+1):(-n_points-1)] -
                      xs[:, (n_points):(-n_points-2)])
    for n in range(1, n_points):
        ptv = ptv + np.absolute(xs[:, (n_points+1+n):(-n_points-1+n)] -
                                xs[:, (n_points+2+n):(-n_points+n)])
        ntv = ntv + np.absolute(xs[:, (n_points+1-n):(-n_points-1-n)] -
                                xs[:, (n_points-n):(-n_points-2-n)])

    if axis:
        return ptv, ntv
    else:
        return ptv.T, ntv.T


def _gibbs_removal_1d(x, axis=0, n_points=3):
    """ Suppresses Gibbs ringing along a given axis using fourier sub-shifts.

    Parameters
    ----------
    x : 2D ndarray
        Matrix x.
    axis : int (0 or 1)
        Axis in which Gibbs oscillations will be suppressed. Default is set
        to 0
    n_points : int, optional
        Number of neighbours to access local TV (see note). Default is set to
        3.

    Returns
    -------
    xc : 2D ndarray
        Matrix with suppressed Gibbs oscillations along the given axis.

    Note
    ----
    This function suppresses the effects of Gibbs oscillations based on the
    analysis of local total variation (TV). Although artefact correction is
    done based on two adjanced points for each voxel, total variation should be
    accessed in a larger range of neighbours. The number of the neighbours to
    be considered in TV calculation can be adjusted using parameter n_points.
    """
    ssamp = np.linspace(0.02, 0.9, num=45)

    xs = x.copy() if axis else x.T.copy()

    # TV for shift zero (baseline)
    tvr, tvl = _image_tv(xs, axis=1, n_points=n_points)
    tvp = np.minimum(tvr, tvl)
    tvn = tvp.copy()

    # Find optimal shift for gibbs removal
    isp = xs.copy()
    isn = xs.copy()
    sp = np.zeros(xs.shape)
    sn = np.zeros(xs.shape)
    N = xs.shape[1]
    c = np.fft.fftshift(np.fft.fft2(xs))
    k = np.linspace(-N/2, N/2-1, num=N)
    k = (2.0j * np.pi * k) / N
    for s in ssamp:
        # Access positive shift for given s
        img_p = abs(np.fft.ifft2(np.fft.fftshift(c * np.exp(k*s))))
        tvsr, tvsl = _image_tv(img_p, axis=1, n_points=n_points)
        tvs_p = np.minimum(tvsr, tvsl)

        # Access negative shift for given s
        img_n = abs(np.fft.ifft2(np.fft.fftshift(c * np.exp(-k*s))))
        tvsr, tvsl = _image_tv(img_n, axis=1, n_points=n_points)
        tvs_n = np.minimum(tvsr, tvsl)

        # Update positive shift params
        isp[tvp > tvs_p] = img_p[tvp > tvs_p]
        sp[tvp > tvs_p] = s
        tvp[tvp > tvs_p] = tvs_p[tvp > tvs_p]

        # Update negative shift params
        isn[tvn > tvs_n] = img_n[tvn > tvs_n]
        sn[tvn > tvs_n] = s
        tvn[tvn > tvs_n] = tvs_n[tvn > tvs_n]

    # check non-zero sub-voxel shifts
    idx = np.nonzero(sp + sn)

    # use positive and negative optimal sub-voxel shifts to interpolate to
    # original grid points
    xs[idx] = (isp[idx] - isn[idx])/(sp[idx] + sn[idx])*sn[idx] + isn[idx]

    return xs if axis else xs.T


def _weights(shape):
    """ Computes the weights necessary to combine two images processed by
    the 1D Gibbs removal procedure along two different axes [1]_.

    Parameters
    ----------
    shape : tuple
        shape of the image

    Returns
    -------
    G0 : 2D ndarray
        Weights for the image corrected along axis 0.
    G1 : 2D ndarray
        Weights for the image corrected along axis 1.

    References
    ----------
    .. [1] <NAME>, <NAME>, <NAME>, <NAME>. Gibbs-ringing artifact
           removal based on local subvoxel-shifts. Magn Reson Med. 2016
           doi: 10.1002/mrm.26054.
    """
    G0 = np.zeros(shape)
    G1 = np.zeros(shape)
    k0 = np.linspace(-np.pi, np.pi, num=shape[0])
    k1 = np.linspace(-np.pi, np.pi, num=shape[1])

    # Middle points
    K1, K0 = np.meshgrid(k1[1:-1], k0[1:-1])
    cosk0 = 1.0 + np.cos(K0)
    cosk1 = 1.0 + np.cos(K1)
    G1[1:-1, 1:-1] = cosk0 / (cosk0 + cosk1)
    G0[1:-1, 1:-1] = cosk1 / (cosk0 + cosk1)

    # Boundaries
    G1[1:-1, 0] = G1[1:-1, -1] = 1
    G1[0, 0] = G1[-1, -1] = G1[0, -1] = G1[-1, 0] = 1/2
    G0[0, 1:-1] = G0[-1, 1:-1] = 1
    G0[0, 0] = G0[-1, -1] = G0[0, -1] = G0[-1, 0] = 1/2

    return G0, G1


def _gibbs_removal_2d(image, n_points=3, G0=None, G1=None):
{'fraction_non_alphanumeric': 0.13624029194348275, 'fraction_numerical': 0.03144006737157294, 'mean_word_length': 2.1078801977318986, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 31, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1510107', 'n_tokens_mistral': 4123, 'n_tokens_neox': 3899, 'n_words': 1146}