def project(obj, theta, center=None, emission=True, sinogram_order=False, ncore=None, nchunk=None):
    """
    Project x-rays through a given 3D object.

    Parameters
    ----------
    obj : ndarray
        Voxelized 3D object.
    theta : array
        Projection angles in radian.
    center: array, optional
        Location of rotation axis.
    emission : bool, optional
        Determines whether output data is emission or transmission type.
    sinogram_order: bool, optional
        Determins whether output data is a stack of sinograms (True, y-axis first axis) 
        or a stack of radiographs (False, theta first axis).
    ncore : int, optional
        Number of cores that will be assigned to jobs.
    nchunk : int, optional
        Chunk size for each core.

    Returns
    -------
    ndarray
        3D tomographic data.
    """
    obj = dtype.as_float32(obj)
    theta = dtype.as_float32(theta)

    # Estimate data dimensions.
    oy, ox, oz = obj.shape
    dt = theta.size
    dy = oy
    dx = _round_to_even(np.sqrt(ox * ox + oz * oz) + 2)
    shape = dy, dt, dx
    tomo = dtype.empty_shared_array(shape)
    tomo[:] = 0.0
    center = get_center(shape, center)

    tomo = mproc.distribute_jobs(
        (obj, center, tomo),
        func=extern.c_project,
        args=(theta,),
        axis=0,
        ncore=ncore,
        nchunk=nchunk)
    # NOTE: returns sinogram order with emmission=True
    if not emission:
        # convert data to be transmission type
        np.exp(-tomo, tomo)
    if not sinogram_order:
        # rotate to radiograph order
        tomo = np.swapaxes(tomo, 0, 1) #doesn't copy data
        # copy data to sharedmem
        tomo = dtype.as_sharedmem(tomo, copy=True)
        
    return tomo

def find_center(
        tomo, theta, ind=None, emission=True, init=None,
        tol=0.5, mask=True, ratio=1.):
    """
    Find rotation axis location.

    The function exploits systematic artifacts in reconstructed images
    due to shifts in the rotation center. It uses image entropy
    as the error metric and ''Nelder-Mead'' routine (of the scipy
    optimization module) as the optimizer :cite:`Donath:06`.

    Parameters
    ----------
    tomo : ndarray
        3D tomographic data.
    theta : array
        Projection angles in radian.
    ind : int, optional
        Index of the slice to be used for reconstruction.
    emission : bool, optional
        Determines whether data is emission or transmission type.
    init : float
        Initial guess for the center.
    tol : scalar
        Desired sub-pixel accuracy.
    mask : bool, optional
        If ``True``, apply a circular mask to the reconstructed image to
        limit the analysis into a circular region.
    ratio : float, optional
        The ratio of the radius of the circular mask to the edge of the
        reconstructed image.

    Returns
    -------
    float
        Rotation axis location.
    """
    tomo = dtype.as_float32(tomo)
    theta = dtype.as_float32(theta)

    if ind is None:
        ind = tomo.shape[1] // 2
    if init is None:
        init = tomo.shape[2] // 2
    print (ind, tomo.shape)

    hmin, hmax = _adjust_hist_limits(
        tomo[:, ind:ind + 1, :], theta, ind, mask, emission)

    # Magic is ready to happen...
    res = minimize(
        _find_center_cost, init,
        args=(tomo, theta, ind, hmin, hmax, mask, ratio, emission),
        method='Nelder-Mead',
        tol=tol)
    return res.x

def sobel_filter(arr, axis=0, ncore=None):
    """
    Apply Sobel filter to 3D array along specified axis.

    Parameters
    ----------
    arr : ndarray
        Input array.
    axis : int, optional
        Axis along which sobel filtering is performed.
    ncore : int, optional
        Number of cores that will be assigned to jobs.

    Returns
    -------
    ndarray
        3D array of same shape as input.
    """
    arr = dtype.as_float32(arr)
    arr = mproc.distribute_jobs(
        arr,
        func=filters.sobel,
        axis=axis,
        ncore=ncore,
        nchunk=0)
    return arr

def circ_mask(arr, axis, ratio=1, val=0.):
    """
    Apply circular mask to a 3D array.

    Parameters
    ----------
    arr : ndarray
        Arbitrary 3D array.
    axis : int
        Axis along which mask will be performed.
    ratio : int, optional
        Ratio of the mask's diameter in pixels to
        the smallest edge size along given axis.
    val : int, optional
        Value for the masked region.

    Returns
    -------
    ndarray
        Masked array.
    """
    arr = dtype.as_float32(arr)
    _arr = arr.swapaxes(0, axis)
    dx, dy, dz = _arr.shape
    mask = _get_mask(dy, dz, ratio)
    for m in range(dx):
        _arr[m, ~mask] = val
    return _arr.swapaxes(0, axis)

def normalize_roi(tomo, roi=[0, 0, 10, 10], ncore=None, nchunk=None):
    """
    Normalize raw projection data using an average of a selected window
    on projection images.

    Parameters
    ----------
    tomo : ndarray
        3D tomographic data.
    roi: list of int, optional
        [top-left, top-right, bottom-left, bottom-right] pixel coordinates.
    ncore : int, optional
        Number of cores that will be assigned to jobs.
    nchunk : int, optional
        Chunk size for each core.

    Returns
    -------
    ndarray
        Normalized 3D tomographic data.
    """
    tomo = dtype.as_float32(tomo)

    arr = mproc.distribute_jobs(
        tomo,
        func=_normalize_roi,
        args=(roi, ),
        axis=0,
        ncore=ncore,
        nchunk=nchunk)
    return arr

def median_filter(arr, size=3, axis=0, ncore=None):
    """
    Apply median filter to 3D array along specified axis.

    Parameters
    ----------
    arr : ndarray
        Input array.
    size : int, optional
        The size of the filter.
    axis : int, optional
        Axis along which median filtering is performed.
    ncore : int, optional
        Number of cores that will be assigned to jobs.

    Returns
    -------
    ndarray
        Median filtered 3D array.
    """
    arr = dtype.as_float32(arr)
    arr = mproc.distribute_jobs(
        arr,
        func=filters.median_filter,
        args=((size, size),),
        axis=axis,
        ncore=ncore,
        nchunk=0)
    return arr

def median_filter(arr, size=3, axis=0, ncore=None):
    """
    Apply median filter to 3D array along specified axis.

    Parameters
    ----------
    arr : ndarray
        Input array.
    size : int, optional
        The size of the filter.
    axis : int, optional
        Axis along which median filtering is performed.
    ncore : int, optional
        Number of cores that will be assigned to jobs.

    Returns
    -------
    ndarray
        Median filtered 3D array.
    """
    arr = dtype.as_float32(arr)
    out = np.empty_like(arr)

    if ncore is None:
        ncore = mproc.mp.cpu_count()

    with cf.ThreadPoolExecutor(ncore) as e:
        slc = [slice(None)]*arr.ndim
        for i in range(arr.shape[axis]):
            slc[axis] = i
            e.submit(filters.median_filter, arr[tuple(slc)], size=(size, size),
                     output=out[tuple(slc)])
    return out

def normalize_bg(tomo, air=1, ncore=None, nchunk=None):
    """
    Normalize 3D tomgraphy data based on background intensity.

    Weight sinogram such that the left and right image boundaries
    (i.e., typically the air region around the object) are set to one
    and all intermediate values are scaled linearly.

    Parameters
    ----------
    tomo : ndarray
        3D tomographic data.
    air : int, optional
        Number of pixels at each boundary to calculate the scaling factor.
    ncore : int, optional
        Number of cores that will be assigned to jobs.
    nchunk : int, optional
        Chunk size for each core.

    Returns
    -------
    ndarray
        Corrected 3D tomographic data.
    """
    tomo = dtype.as_float32(tomo)
    air = dtype.as_int32(air)

    arr = mproc.distribute_jobs(
        tomo,
        func=extern.c_normalize_bg,
        args=(air,),
        axis=0,
        ncore=ncore,
        nchunk=nchunk)
    return arr

def remove_nan(arr, val=0., ncore=None):
    """
    Replace NaN values in array with a given value.

    Parameters
    ----------
    arr : ndarray
        Input array.
    val : float, optional
        Values to be replaced with NaN values in array.
    ncore : int, optional
        Number of cores that will be assigned to jobs.

    Returns
    -------
    ndarray
       Corrected array.
    """
    arr = dtype.as_float32(arr)
    val = np.float32(val)

    with mproc.set_numexpr_threads(ncore):
        ne.evaluate('where(arr!=arr, val, arr)', out=arr)

    return arr

def sobel_filter(arr, axis=0, ncore=None):
    """
    Apply Sobel filter to 3D array along specified axis.

    Parameters
    ----------
    arr : ndarray
        Input array.
    axis : int, optional
        Axis along which sobel filtering is performed.
    ncore : int, optional
        Number of cores that will be assigned to jobs.

    Returns
    -------
    ndarray
        3D array of same shape as input.
    """
    arr = dtype.as_float32(arr)
    out = np.empty_like(arr)

    if ncore is None:
        ncore = mproc.mp.cpu_count()

    with cf.ThreadPoolExecutor(ncore) as e:
        slc = [slice(None)]*arr.ndim
        for i in range(arr.shape[axis]):
            slc[axis] = i
            e.submit(filters.sobel, arr[slc], output=out[slc])
    return out

def circ_mask(arr, axis, ratio=1, val=0., ncore=None):
    """
    Apply circular mask to a 3D array.

    Parameters
    ----------
    arr : ndarray
            Arbitrary 3D array.
    axis : int
        Axis along which mask will be performed.
    ratio : int, optional
        Ratio of the mask's diameter in pixels to
        the smallest edge size along given axis.
    val : int, optional
        Value for the masked region.

    Returns
    -------
    ndarray
        Masked array.
    """
    arr = dtype.as_float32(arr)
    val = np.float32(val)
    _arr = arr.swapaxes(0, axis)
    dx, dy, dz = _arr.shape
    mask = _get_mask(dy, dz, ratio)

    with mproc.set_numexpr_threads(ncore):
        ne.evaluate('where(mask, _arr, val)', out=_arr)

    return _arr.swapaxes(0, axis)

def vector(tomo, theta, center=None, num_iter=1):
    tomo = dtype.as_float32(tomo)
    theta = dtype.as_float32(theta)

    # Initialize tomography data.
    tomo = init_tomo(tomo, sinogram_order=False, sharedmem=False)

    recon_shape = (tomo.shape[0], tomo.shape[2], tomo.shape[2])
    recon1 = np.zeros(recon_shape, dtype=np.float32)
    recon2 = np.zeros(recon_shape, dtype=np.float32)

    center_arr = get_center(tomo.shape, center)

    extern.c_vector(tomo, center_arr, recon1, recon2, theta, 
        num_gridx=tomo.shape[2], num_gridy=tomo.shape[2], num_iter=num_iter)
    return recon1, recon2

def remove_stripe_sf(tomo, size=5, ncore=None, nchunk=None):
    """
    Normalize raw projection data using a smoothing filter approach.

    Parameters
    ----------
    tomo : ndarray
        3D tomographic data.
    size : int, optional
        Size of the smoothing filter.
    ncore : int, optional
        Number of cores that will be assigned to jobs.
    nchunk : int, optional
        Chunk size for each core.

    Returns
    -------
    ndarray
        Corrected 3D tomographic data.
    """
    tomo = dtype.as_float32(tomo)
    dx, dy, dz = tomo.shape
    arr = mproc.distribute_jobs(
        tomo,
        func=extern.c_remove_stripe_sf,
        args=(dx, dy, dz, size),
        axis=1,
        ncore=ncore,
        nchunk=nchunk)
    return arr

def remove_outlier(arr, dif, size=3, axis=0, ncore=None):
    """
    Remove high intensity bright spots from a 3D array along specified
    dimension.

    Parameters
    ----------
    arr : ndarray
        Input array.
    dif : float
        Expected difference value between outlier value and
        the median value of the array.
    size : int
        Size of the median filter.
    axis : int, optional
        Axis along which median filtering is performed.
    ncore : int, optional
        Number of cores that will be assigned to jobs.

    Returns
    -------
    ndarray
       Corrected array.
    """
    arr = dtype.as_float32(arr)
    arr = mproc.distribute_jobs(
        arr,
        func=_remove_outlier_from_img,
        args=(dif, size),
        axis=axis,
        ncore=ncore,
        nchunk=0)
    return arr

def find_center_vo(tomo, ind=None, smin=-50, smax=50, srad=6, step=0.25,
                   ratio=0.5, drop=20):
    """
    Find rotation axis location using Nghia Vo's method. :cite:`Vo:14`.

    Parameters
    ----------
    tomo : ndarray
        3D tomographic data.
    ind : int, optional
        Index of the slice to be used for reconstruction.
    smin, smax : int, optional
        Coarse search radius. Reference to the horizontal center of the sinogram.
    srad : float, optional
        Fine search radius.
    step : float, optional
        Step of fine searching.
    ratio : float, optional
        The ratio between the FOV of the camera and the size of object.
        It's used to generate the mask.
    drop : int, optional
        Drop lines around vertical center of the mask.

    Returns
    -------
    float
        Rotation axis location.
    """
    tomo = dtype.as_float32(tomo)
    (depth, height, width) = tomo.shape
    if ind is None:
        ind = height // 2
        if height > 10:
            # Averaging sinograms to improve SNR
            _tomo = np.mean(tomo[:, ind - 5:ind + 5, :], axis=1)
        else:
            _tomo = tomo[:, ind, :]
    else:
        _tomo = tomo[:, ind, :]

    # Denoising
    # There's a critical reason to use different window sizes
    # between coarse and fine search.
    _tomo_cs = ndimage.filters.gaussian_filter(_tomo, (3, 1))
    _tomo_fs = ndimage.filters.gaussian_filter(_tomo, (2, 2))

    # Coarse and fine searches for finding the rotation center.
    if _tomo.shape[0] * _tomo.shape[1] > 4e6:  # If data is large (>2kx2k)
        _tomo_coarse = downsample(
            np.expand_dims(_tomo_cs, 1), level=2)[:, 0, :]
        init_cen = _search_coarse(
            _tomo_coarse, smin / 4.0, smax / 4.0, ratio, drop)
        fine_cen = _search_fine(_tomo_fs, srad, step,
                                init_cen * 4, ratio, drop)
    else:
        init_cen = _search_coarse(_tomo_cs, smin, smax, ratio, drop)
        fine_cen = _search_fine(_tomo_fs, srad, step, init_cen, ratio, drop)

    logger.debug('Rotation center search finished: %i', fine_cen)
    return fine_cen

def remove_ring(rec, center_x=None, center_y=None, thresh=300.0,
                thresh_max=300.0, thresh_min=-100.0, theta_min=30,
                rwidth=30, ncore=None, nchunk=None):
    """
    Remove ring artifacts from images in the reconstructed domain.
    Descriptions of parameters need to be more clear for sure.

    Parameters
    ----------
    arr : ndarray
        Array of reconstruction data
    center_x : float, optional
        abscissa location of center of rotation
    center_y : float, optional
        ordinate location of center of rotation
    thresh : float, optional
        maximum value of an offset due to a ring artifact
    thresh_max : float, optional
        max value for portion of image to filter
    thresh_min : float, optional
        min value for portion of image to filer
    theta_min : int, optional
        minimum angle in degrees (int) to be considered ring artifact
    rwidth : int, optional
        Maximum width of the rings to be filtered in pixels
    ncore : int, optional
        Number of cores that will be assigned to jobs.
    nchunk : int, optional
        Chunk size for each core.

    Returns
    -------
    ndarray
        Corrected reconstruction data
    """

    rec = dtype.as_float32(rec)

    dz, dy, dx = rec.shape

    if center_x is None:
        center_x = (dx - 1.0)/2.0
    if center_y is None:
        center_y = (dy - 1.0)/2.0

    args = (center_x, center_y, dx, dy, dz, thresh_max, thresh_min,
            thresh, theta_min, rwidth)

    rec = mproc.distribute_jobs(
        rec,
        func=extern.c_remove_ring,
        args=args,
        axis=0,
        ncore=ncore,
        nchunk=nchunk)
    return rec

def project(obj, theta, center=None, ncore=None, nchunk=None):
    """
    Project x-rays through a given 3D object.

    Parameters
    ----------
    obj : ndarray
        Voxelized 3D object.
    theta : array
        Projection angles in radian.
    center: array, optional
        Location of rotation axis.
    ncore : int, optional
        Number of cores that will be assigned to jobs.
    nchunk : int, optional
        Chunk size for each core.

    Returns
    -------
    ndarray
        3D tomographic data.
    """
    obj = dtype.as_float32(obj)
    theta = dtype.as_float32(theta)

    # Estimate data dimensions.
    ox, oy, oz = obj.shape
    dx = theta.size
    dy = ox
    dz = np.ceil(np.sqrt(oy * oy + oz * oz)).astype('int')
    shape = dx, dy, dz
    tomo = np.zeros(shape, dtype='float32')
    center = get_center(shape, center)

    mproc.init_obj(obj)
    arr = mproc.distribute_jobs(
        tomo,
        func=extern.c_project,
        args=(ox, oy, oz, theta, center, dx, dy, dz),
        axis=0,
        ncore=ncore,
        nchunk=nchunk)
    return arr

def _sample(arr, level, axis, mode):
    arr = dtype.as_float32(arr.copy())
    dx, dy, dz = arr.shape
    # Determine the new size, dim, of the down-/up-sampled dimension
    if mode == 0:
        dim = int(arr.shape[axis] / np.power(2, level))
    if mode == 1:
        dim = int(arr.shape[axis] * np.power(2, level))

    out = _init_out(arr, axis, dim)
    return extern.c_sample(mode, arr, dx, dy, dz, level, axis, out)

def init_tomo(tomo, sinogram_order, sharedmem=True):
    tomo = dtype.as_float32(tomo)
    if not sinogram_order:
        tomo = np.swapaxes(tomo, 0, 1) #doesn't copy data
    if sharedmem:
        # copy data to sharedmem (if not already or not contiguous)
        tomo = dtype.as_sharedmem(tomo, copy=not dtype.is_contiguous(tomo))
    else:
        # ensure contiguous
        tomo = np.require(tomo, requirements="AC")
    return tomo

def _sample(arr, level, axis, mode):
    arr = dtype.as_float32(arr)
    dx, dy, dz = arr.shape

    if mode == 0:
        dim = arr.shape[axis] / np.power(2, level)
    if mode == 1:
        dim = arr.shape[axis] * np.power(2, level)

    out = _init_out(arr, axis, dim)
    return extern.c_sample(mode, arr, dx, dy, dz, level, axis, out)