# Path: pyem/vop/vop_numba.py
# @numba.jit(cache=True, nopython=True, nogil=True)
# def fill_ft(ft, ftc, rmax, normfft=1):
#     rmax2 = rmax ** 2
#     for k in range(ft.shape[0]):
#         kp = k if k < ft.shape[2] else k - ft.shape[0]
#         for i in range(ft.shape[1]):
#             ip = i if i < ft.shape[2] else i - ft.shape[1]
#             for j in range(ft.shape[2]):
#                 jp = j
#                 r2 = ip**2 + jp**2 + kp**2
#                 if r2 <= rmax2:
#                     ftc[kp + ftc.shape[0]//2, ip + ftc.shape[1]//2, jp] = ft[k, i, j] * normfft


import numpy as np
import numpy.ma as ma
from scipy.ndimage import map_coordinates
from pyfftw.interfaces.numpy_fft import rfftn
from .vop_numba import fill_ft
def interpolate_slice(f3d, rot, pfac=2, size=None):
    nhalf = f3d.shape[0] / 2
    if size is None:
        phalf = nhalf
    else:
        phalf = size / 2
    qot = rot * pfac  # Scaling!
    px, py, pz = np.meshgrid(np.arange(-phalf, phalf), np.arange(-phalf, phalf), 0)
    pr = np.sqrt(px ** 2 + py ** 2 + pz ** 2)
    pcoords = np.vstack([px.reshape(-1), py.reshape(-1), pz.reshape(-1)])
    mcoords = qot.T.dot(pcoords)
    mcoords = mcoords[:, pr.reshape(-1) < nhalf]
    pvals = map_coordinates(np.real(f3d), mcoords, order=1, mode="wrap") + \
             1j * map_coordinates(np.imag(f3d), mcoords, order=1, mode="wrap")
    pslice = np.zeros(pr.shape, dtype=np.complex)
    pslice[pr < nhalf] = pvals
    return pslice


def vol_ft(vol, pfac=2, threads=1, normfft=1):
    """ Returns a centered, Nyquist-limited, zero-padded, interpolation-ready 3D Fourier transform.
    :param vol: Volume to be Fourier transformed.
    :param pfac: Size factor for zero-padding.
    :param threads: Number of threads for pyFFTW.
    :param normfft: Normalization constant for Fourier transform.
    """
    vol = grid_correct(vol, pfac=pfac, order=1)
    padvol = np.pad(vol, int((vol.shape[0] * pfac - vol.shape[0]) // 2), "constant")
    ft = rfftn(np.fft.ifftshift(padvol), padvol.shape, threads=threads)
    ftc = np.zeros((ft.shape[0] + 3, ft.shape[1] + 3, ft.shape[2]), dtype=ft.dtype)

# Path: pyem/vop/vop_numba.py
# @numba.jit(cache=True, nopython=True, nogil=True)
# def fill_ft(ft, ftc, rmax, normfft=1):
#     rmax2 = rmax ** 2
#     for k in range(ft.shape[0]):
#         kp = k if k < ft.shape[2] else k - ft.shape[0]
#         for i in range(ft.shape[1]):
#             ip = i if i < ft.shape[2] else i - ft.shape[1]
#             for j in range(ft.shape[2]):
#                 jp = j
#                 r2 = ip**2 + jp**2 + kp**2
#                 if r2 <= rmax2:
#                     ftc[kp + ftc.shape[0]//2, ip + ftc.shape[1]//2, jp] = ft[k, i, j] * normfft


# Path: pyem/vop/vop.py
import numpy as np
import numpy.ma as ma
from scipy.ndimage import map_coordinates
from pyfftw.interfaces.numpy_fft import rfftn
from .vop_numba import fill_ft

def interpolate_slice(f3d, rot, pfac=2, size=None):
    nhalf = f3d.shape[0] / 2
    if size is None:
        phalf = nhalf
    else:
        phalf = size / 2
    qot = rot * pfac  # Scaling!
    px, py, pz = np.meshgrid(np.arange(-phalf, phalf), np.arange(-phalf, phalf), 0)
    pr = np.sqrt(px ** 2 + py ** 2 + pz ** 2)
    pcoords = np.vstack([px.reshape(-1), py.reshape(-1), pz.reshape(-1)])
    mcoords = qot.T.dot(pcoords)
    mcoords = mcoords[:, pr.reshape(-1) < nhalf]
    pvals = map_coordinates(np.real(f3d), mcoords, order=1, mode="wrap") + \
             1j * map_coordinates(np.imag(f3d), mcoords, order=1, mode="wrap")
    pslice = np.zeros(pr.shape, dtype=np.complex)
    pslice[pr < nhalf] = pvals
    return pslice


def vol_ft(vol, pfac=2, threads=1, normfft=1):
    """ Returns a centered, Nyquist-limited, zero-padded, interpolation-ready 3D Fourier transform.
    :param vol: Volume to be Fourier transformed.
    :param pfac: Size factor for zero-padding.
    :param threads: Number of threads for pyFFTW.
    :param normfft: Normalization constant for Fourier transform.
    """
    vol = grid_correct(vol, pfac=pfac, order=1)
    padvol = np.pad(vol, int((vol.shape[0] * pfac - vol.shape[0]) // 2), "constant")
    ft = rfftn(np.fft.ifftshift(padvol), padvol.shape, threads=threads)
    ftc = np.zeros((ft.shape[0] + 3, ft.shape[1] + 3, ft.shape[2]), dtype=ft.dtype)

    fill_ft(ft, ftc, vol.shape[0], normfft=normfft)