def test_gridrec(self):
     assert_allclose(
         recon(self.prj, self.ang, algorithm='gridrec', filter_name='none'),
         read_file('gridrec_none.npy'), rtol=1e-2)
     assert_allclose(
         recon(self.prj, self.ang, algorithm='gridrec', filter_name='shepp'),
         read_file('gridrec_shepp.npy'), rtol=1e-2)
     assert_allclose(
         recon(self.prj, self.ang, algorithm='gridrec', filter_name='cosine'),
         read_file('gridrec_cosine.npy'), rtol=1e-2)
     assert_allclose(
         recon(self.prj, self.ang, algorithm='gridrec', filter_name='hann'),
         read_file('gridrec_hann.npy'), rtol=1e-2)
     assert_allclose(
         recon(self.prj, self.ang, algorithm='gridrec', filter_name='hamming'),
         read_file('gridrec_hamming.npy'), rtol=1e-2)
     assert_allclose(
         recon(self.prj, self.ang, algorithm='gridrec', filter_name='ramlak'),
         read_file('gridrec_ramlak.npy'), rtol=1e-2)
     assert_allclose(
         recon(self.prj, self.ang, algorithm='gridrec', filter_name='parzen'),
         read_file('gridrec_parzen.npy'), rtol=1e-2)
     assert_allclose(
         recon(self.prj, self.ang, algorithm='gridrec', filter_name='butterworth'),
         read_file('gridrec_butterworth.npy'), rtol=1e-2)

def _adjust_hist_limits(tomo, theta, ind, mask, emission):
    # Make an initial reconstruction to adjust histogram limits.
    rec = recon(tomo, theta, emission=emission, algorithm='gridrec')

    # Apply circular mask.
    if mask is True:
        rec = circ_mask(rec, axis=0)

    # Adjust histogram boundaries according to reconstruction.
    return _adjust_hist_min(rec.min()), _adjust_hist_max(rec.max())

def _find_center_cost(
        center, tomo, theta, ind, hmin, hmax, mask, ratio, emission):
    """
    Cost function used for the ``find_center`` routine.
    """
    logger.info('trying center: %s', center)
    center = np.array(center, dtype='float32')
    rec = recon(
        tomo[:, ind:ind + 1, :], theta, center, emission=emission, algorithm='gridrec')

    if mask is True:
        rec = circ_mask(rec, axis=0)

    hist, e = np.histogram(rec, bins=64, range=[hmin, hmax])
    hist = hist.astype('float32') / rec.size + 1e-12
    return -np.dot(hist, np.log2(hist))

def _find_center_cost(
        center, tomo_ind, theta, hmin, hmax, mask, ratio, 
        sinogram_order=False):
    """
    Cost function used for the ``find_center`` routine.
    """
    logger.info('Trying rotation center: %s', center)
    center = np.array(center, dtype='float32')
    rec = recon(
        tomo_ind, theta, center,
        sinogram_order=sinogram_order, algorithm='gridrec')

    if mask is True:
        rec = circ_mask(rec, axis=0)

    hist, e = np.histogram(rec, bins=64, range=[hmin, hmax])
    hist = hist.astype('float32') / rec.size + 1e-12
    val = -np.dot(hist, np.log2(hist))
    logger.info("Function value = %f"%val)    
    return val

 def test_gridrec_custom(self):
     assert_allclose(
         recon(self.prj, self.ang, algorithm='gridrec', filter_name='none'),
         recon(
             self.prj, self.ang, algorithm='gridrec', filter_name='custom',
             filter_par=np.ones(self.prj.shape[-1], dtype=np.float32)))

 def test_fbp(self):
     assert_allclose(
         recon(self.prj, self.ang, algorithm='fbp'),
         read_file('fbp.npy'), rtol=1e-2)

 def test_sirt(self):
     # FIXME: Make separate tests for each back-end
     os.environ["TOMOPY_USE_C_SIRT"] = "1"
     r_sirt = recon(self.prj, self.ang, algorithm='sirt', num_iter=4)
     c_sirt = read_file('sirt.npy')
     assert_allclose(r_sirt, c_sirt, rtol=1e-2)

 def test_ospml_hybrid(self):
     assert_allclose(
         recon(self.prj, self.ang, algorithm='ospml_hybrid', num_iter=4),
         read_file('ospml_hybrid.npy'), rtol=1e-2)

def write_center(
        tomo, theta, dpath='tmp/center', cen_range=None, ind=None,
        mask=False, ratio=1., sinogram_order=False, algorithm='gridrec', filter_name='parzen'):
    """
    Save images reconstructed with a range of rotation centers.

    Helps finding the rotation center manually by visual inspection of
    images reconstructed with a set of different centers.The output
    images are put into a specified folder and are named by the
    center position corresponding to the image.

    Parameters
    ----------
    tomo : ndarray
        3D tomographic data.
    theta : array
        Projection angles in radian.
    dpath : str, optional
        Folder name to save output images.
    cen_range : list, optional
        [start, end, step] Range of center values.
    ind : int, optional
        Index of the slice to be used for reconstruction.
    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.
    sinogram_order: bool, optional
        Determins whether data is a stack of sinograms (True, y-axis first axis)
        or a stack of radiographs (False, theta first axis).
    algorithm : {str, function}
        One of the following string values.

        'art'
            Algebraic reconstruction technique :cite:`Kak:98`.
        'bart'
            Block algebraic reconstruction technique.
        'fbp'
            Filtered back-projection algorithm.
        'gridrec'
            Fourier grid reconstruction algorithm :cite:`Dowd:99`,
            :cite:`Rivers:06`.
        'mlem'
            Maximum-likelihood expectation maximization algorithm
            :cite:`Dempster:77`.
        'osem'
            Ordered-subset expectation maximization algorithm
            :cite:`Hudson:94`.
        'ospml_hybrid'
            Ordered-subset penalized maximum likelihood algorithm with
            weighted linear and quadratic penalties.
        'ospml_quad'
            Ordered-subset penalized maximum likelihood algorithm with
            quadratic penalties.
        'pml_hybrid'
            Penalized maximum likelihood algorithm with weighted linear
            and quadratic penalties :cite:`Chang:04`.
        'pml_quad'
            Penalized maximum likelihood algorithm with quadratic penalty.
        'sirt'
            Simultaneous algebraic reconstruction technique.
        'tv'
            Total Variation reconstruction technique
            :cite:`Chambolle:11`.
        'grad'
            Gradient descent method with a constant step size

    filter_name : str, optional
        Name of the filter for analytic reconstruction.

        'none'
            No filter.
        'shepp'
            Shepp-Logan filter (default).
        'cosine'
            Cosine filter.
        'hann'
            Cosine filter.
        'hamming'
            Hamming filter.
        'ramlak'
            Ram-Lak filter.
        'parzen'
            Parzen filter.
        'butterworth'
            Butterworth filter.
        'custom'
            A numpy array of size `next_power_of_2(num_detector_columns)/2`
            specifying a custom filter in Fourier domain. The first element
            of the filter should be the zero-frequency component.
        'custom2d'
            A numpy array of size `num_projections*next_power_of_2(num_detector_columns)/2`
            specifying a custom angle-dependent filter in Fourier domain. The first element
            of each filter should be the zero-frequency component.
    """
    tomo = dtype.as_float32(tomo)
    theta = dtype.as_float32(theta)

#.........这里部分代码省略.........

 def test_mlem(self):
     assert_allclose(
         recon(self.prj, self.ang, algorithm='mlem', num_iter=4),
         read_file('mlem.npy'), rtol=1e-2)

def align_seq(
        prj, ang, fdir='.', iters=10, pad=(0, 0),
        blur=True, save=False, debug=True):
    """
    Aligns the projection image stack using the sequential
    re-projection algorithm :cite:`Gursoy:17`.

    Parameters
    ----------
    prj : ndarray
        3D stack of projection images. The first dimension
        is projection axis, second and third dimensions are
        the x- and y-axes of the projection image, respectively.
    ang : ndarray
        Projection angles in radians as an array.
    iters : scalar, optional
        Number of iterations of the algorithm.
    pad : list-like, optional
        Padding for projection images in x and y-axes.
    blur : bool, optional
        Blurs the edge of the image before registration.
    save : bool, optional
        Saves projections and corresponding reconstruction
        for each algorithm iteration.
    debug : book, optional
        Provides debugging info such as iterations and error.

    Returns
    -------
    ndarray
        3D stack of projection images with jitter.
    ndarray
        Error array for each iteration.
    """

    # Needs scaling for skimage float operations.
    prj, scl = scale(prj)

    # Shift arrays
    sx = np.zeros((prj.shape[0]))
    sy = np.zeros((prj.shape[0]))

    conv = np.zeros((iters))

    # Pad images.
    npad = ((0, 0), (pad[1], pad[1]), (pad[0], pad[0]))
    prj = np.pad(prj, npad, mode='constant', constant_values=0)

    # Register each image frame-by-frame.
    for n in range(iters):
        # Reconstruct image.
        rec = recon(prj, ang, algorithm='sirt')

        # Re-project data and obtain simulated data.
        sim = project(rec, ang, pad=False)

        # Blur edges.
        if blur:
            _prj = blur_edges(prj, 0.1, 0.5)
            _sim = blur_edges(sim, 0.1, 0.5)
        else:
            _prj = prj
            _sim = sim

        # Initialize error matrix per iteration.
        err = np.zeros((prj.shape[0]))

        # For each projection
        for m in range(prj.shape[0]):

            # Register current projection in sub-pixel precision
            shift, error, diffphase = register_translation(_prj[m], _sim[m], 2)
            err[m] = np.sqrt(shift[0]*shift[0] + shift[1]*shift[1])
            sx[m] += shift[0]
            sy[m] += shift[1]

            # Register current image with the simulated one
            tform = tf.SimilarityTransform(translation=(shift[1], shift[0]))
            prj[m] = tf.warp(prj[m], tform, order=5)

        if debug:
            print('iter=' + str(n) + ', err=' + str(np.linalg.norm(err)))
            conv[n] = np.linalg.norm(err)

        if save:
            dxchange.write_tiff(prj, fdir + '/tmp/iters/prj/prj')
            dxchange.write_tiff(sim, fdir + '/tmp/iters/sim/sim')
            dxchange.write_tiff(rec, fdir + '/tmp/iters/rec/rec')

    # Re-normalize data
    prj *= scl
    return prj, sx, sy, conv

 def test_art(self):
     os.environ["TOMOPY_USE_C_ART"] = "1"
     assert_allclose(
         recon(self.prj, self.ang, algorithm='art', num_iter=4),
         read_file('art.npy'), rtol=1e-2)

 def test_tv(self):
     assert_allclose(
         recon(self.prj, self.ang, algorithm='tv', num_iter=4),
         read_file('tv.npy'), rtol=1e-2)

def write_center(
        tomo, theta, dpath='tmp/center', cen_range=None, ind=None,
        mask=False, ratio=1., sinogram_order=False):
    """
    Save images reconstructed with a range of rotation centers.

    Helps finding the rotation center manually by visual inspection of
    images reconstructed with a set of different centers.The output
    images are put into a specified folder and are named by the
    center position corresponding to the image.

    Parameters
    ----------
    tomo : ndarray
        3D tomographic data.
    theta : array
        Projection angles in radian.
    dpath : str, optional
        Folder name to save output images.
    cen_range : list, optional
        [start, end, step] Range of center values.
    ind : int, optional
        Index of the slice to be used for reconstruction.
    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.
    sinogram_order: bool, optional
        Determins whether data is a stack of sinograms (True, y-axis first axis) 
        or a stack of radiographs (False, theta first axis).        
    """
    tomo = dtype.as_float32(tomo)
    theta = dtype.as_float32(theta)

    if sinogram_order:
        dy, dt, dx = tomo.shape
    else:
        dt, dy, dx = tomo.shape
    if ind is None:
        ind = dy // 2
    if cen_range is None:
        center = np.arange(dx / 2 - 5, dx / 2 + 5, 0.5)
    else:
        center = np.arange(*cen_range)

    stack = dtype.empty_shared_array((len(center), dt, dx))
        
    for m in range(center.size):
        if sinogram_order:
            stack[m] = tomo[ind]
        else:
            stack[m] = tomo[:, ind, :]

    # Reconstruct the same slice with a range of centers.
    rec = recon(stack, 
                theta, 
                center=center, 
                sinogram_order=True, 
                algorithm='gridrec',
                nchunk=1)

    # Apply circular mask.
    if mask is True:
        rec = circ_mask(rec, axis=0)

    # Save images to a temporary folder.
    for m in range(len(center)):
        fname = os.path.join(
            dpath, str('{0:.2f}'.format(center[m]) + '.tiff'))
        dxchange.write_tiff(rec[m], fname=fname, overwrite=True)

 def test_pml_quad(self):
     assert_allclose(
         recon(self.prj, self.ang, algorithm='pml_quad', num_iter=4),
         read_file('pml_quad.npy'), rtol=1e-2)

def align_seq(
        prj, ang, fdir='.', iters=10, pad=(0, 0),
        blur=True, center=None, algorithm='sirt',
        upsample_factor=10, rin=0.5, rout=0.8,
        save=False, debug=True):
    """
    Aligns the projection image stack using the sequential
    re-projection algorithm :cite:`Gursoy:17`.

    Parameters
    ----------
    prj : ndarray
        3D stack of projection images. The first dimension
        is projection axis, second and third dimensions are
        the x- and y-axes of the projection image, respectively.
    ang : ndarray
        Projection angles in radians as an array.
    iters : scalar, optional
        Number of iterations of the algorithm.
    pad : list-like, optional
        Padding for projection images in x and y-axes.
    blur : bool, optional
        Blurs the edge of the image before registration.
    center: array, optional
        Location of rotation axis.
    algorithm : {str, function}
        One of the following string values.

        'art'
            Algebraic reconstruction technique :cite:`Kak:98`.
        'gridrec'
            Fourier grid reconstruction algorithm :cite:`Dowd:99`,
            :cite:`Rivers:06`.
        'mlem'
            Maximum-likelihood expectation maximization algorithm
            :cite:`Dempster:77`.
        'sirt'
            Simultaneous algebraic reconstruction technique.
        'tv'
            Total Variation reconstruction technique
            :cite:`Chambolle:11`.
        'grad'
            Gradient descent method with a constant step size

    upsample_factor : integer, optional
        The upsampling factor. Registration accuracy is
        inversely propotional to upsample_factor. 
    rin : scalar, optional
        The inner radius of blur function. Pixels inside
        rin is set to one.
    rout : scalar, optional
        The outer radius of blur function. Pixels outside
        rout is set to zero.
    save : bool, optional
        Saves projections and corresponding reconstruction
        for each algorithm iteration.
    debug : book, optional
        Provides debugging info such as iterations and error.

    Returns
    -------
    ndarray
        3D stack of projection images with jitter.
    ndarray
        Error array for each iteration.
    """

    # Needs scaling for skimage float operations.
    prj, scl = scale(prj)

    # Shift arrays
    sx = np.zeros((prj.shape[0]))
    sy = np.zeros((prj.shape[0]))

    conv = np.zeros((iters))

    # Pad images.
    npad = ((0, 0), (pad[1], pad[1]), (pad[0], pad[0]))
    prj = np.pad(prj, npad, mode='constant', constant_values=0)

    # Register each image frame-by-frame.
    for n in range(iters):
        # Reconstruct image.
        rec = recon(prj, ang, center=center, algorithm=algorithm)

        # Re-project data and obtain simulated data.
        sim = project(rec, ang, center=center, pad=False)

        # Blur edges.
        if blur:
            _prj = blur_edges(prj, rin, rout)
            _sim = blur_edges(sim, rin, rout)
        else:
            _prj = prj
            _sim = sim

        # Initialize error matrix per iteration.
        err = np.zeros((prj.shape[0]))

        # For each projection
#.........这里部分代码省略.........

 def test_bart(self):
     assert_allclose(
         recon(self.prj, self.ang, algorithm='bart', num_iter=4),
         read_file('bart.npy'), rtol=1e-2)

def write_center(
        tomo, theta, dpath='tmp/center', cen_range=None, ind=None,
        emission=True, mask=False, ratio=1.):
    """
    Save images reconstructed with a range of rotation centers.

    Helps finding the rotation center manually by visual inspection of
    images reconstructed with a set of different centers.The output
    images are put into a specified folder and are named by the
    center position corresponding to the image.

    Parameters
    ----------
    tomo : ndarray
        3D tomographic data.
    theta : array
        Projection angles in radian.
    dpath : str, optional
        Folder name to save output images.
    cen_range : list, optional
        [start, end, step] Range of center values.
    ind : int, optional
        Index of the slice to be used for reconstruction.
    emission : bool, optional
        Determines whether data is emission or transmission type.
    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.
    """
    tomo = dtype.as_float32(tomo)
    theta = dtype.as_float32(theta)

    dx, dy, dz = tomo.shape
    if ind is None:
        ind = dy // 2
    if cen_range is None:
        center = np.arange(dz / 2 - 5, dz / 2 + 5, 0.5)
    if len(cen_range) < 3:
        cen_range[2] = 1
    else:
        center = np.arange(cen_range[0], cen_range[1], cen_range[2] / 2.)

    stack = np.zeros((dx, len(center), dz))
    for m in range(center.size):
        stack[:, m, :] = tomo[:, ind, :]

    # Reconstruct the same slice with a range of centers.
    rec = recon(
        stack, theta, center=center, emission=emission, algorithm='gridrec')

    # Apply circular mask.
    if mask is True:
        rec = circ_mask(rec, axis=0)

    # Save images to a temporary folder.
    for m in range(len(center)):
        if m % 2 == 0:  # 2 slices same bec of gridrec.
            fname = os.path.join(
                dpath, str('{0:.2f}'.format(center[m]) + '.tiff'))
            write_tiff(rec[m:m + 1], fname=fname, overwrite=True)

 def test_mlem(self):
     # FIXME: Make separate tests for each back-end
     os.environ["TOMOPY_USE_C_MLEM"] = "1"
     assert_allclose(
         recon(self.prj, self.ang, algorithm='mlem', num_iter=4),
         read_file('mlem.npy'), rtol=1e-2)