code: import numpy as np
  from skimage.segmentation import random_walker
  from skimage.transform import resize
  from skimage._shared._warnings import expected_warnings
  from skimage._shared import testing
  
  
  # older versions of scipy raise a warning with new NumPy because they use
  # numpy.rank() instead of arr.ndim or numpy.linalg.matrix_rank.
  SCIPY_EXPECTED = 'numpy.linalg.matrix_rank|\A\Z'
  PYAMG_EXPECTED_WARNING = 'pyamg|\A\Z'
  PYAMG_SCIPY_EXPECTED = SCIPY_EXPECTED + '|' + PYAMG_EXPECTED_WARNING
  
  
  def make_2d_syntheticdata(lx, ly=None):
      if ly is None:
          ly = lx
      np.random.seed(1234)
      data = np.zeros((lx, ly)) + 0.1 * np.random.randn(lx, ly)
      small_l = int(lx // 5)
      data[lx // 2 - small_l:lx // 2 + small_l,
           ly // 2 - small_l:ly // 2 + small_l] = 1
      data[lx // 2 - small_l + 1:lx // 2 + small_l - 1,
           ly // 2 - small_l + 1:ly // 2 + small_l - 1] = (
              0.1 * np.random.randn(2 * small_l - 2, 2 * small_l - 2))
      data[lx // 2 - small_l, ly // 2 - small_l // 8:ly // 2 + small_l // 8] = 0
      seeds = np.zeros_like(data)
      seeds[lx // 5, ly // 5] = 1
      seeds[lx // 2 + small_l // 4, ly // 2 - small_l // 4] = 2
      return data, seeds
  
  
  def make_3d_syntheticdata(lx, ly=None, lz=None):
      if ly is None:
          ly = lx
      if lz is None:
          lz = lx
      np.random.seed(1234)
      data = np.zeros((lx, ly, lz)) + 0.1 * np.random.randn(lx, ly, lz)
      small_l = int(lx // 5)
      data[lx // 2 - small_l:lx // 2 + small_l,
           ly // 2 - small_l:ly // 2 + small_l,
           lz // 2 - small_l:lz // 2 + small_l] = 1
      data[lx // 2 - small_l + 1:lx // 2 + small_l - 1,
           ly // 2 - small_l + 1:ly // 2 + small_l - 1,
           lz // 2 - small_l + 1:lz // 2 + small_l - 1] = 0
      # make a hole
      hole_size = np.max([1, small_l // 8])
      data[lx // 2 - small_l,
           ly // 2 - hole_size:ly // 2 + hole_size,
           lz // 2 - hole_size:lz // 2 + hole_size] = 0
      seeds = np.zeros_like(data)
      seeds[lx // 5, ly // 5, lz // 5] = 1
      seeds[lx // 2 + small_l // 4,
            ly // 2 - small_l // 4,
            lz // 2 - small_l // 4] = 2
      return data, seeds
  
  
  def test_2d_bf():
      lx = 70
      ly = 100
      data, labels = make_2d_syntheticdata(lx, ly)
      labels_bf = random_walker(data, labels, beta=90, mode='bf')
      assert (labels_bf[25:45, 40:60] == 2).all()
      assert data.shape == labels.shape
      full_prob_bf = random_walker(data, labels, beta=90, mode='bf',
                                   return_full_prob=True)
      assert (full_prob_bf[1, 25:45, 40:60] >=
              full_prob_bf[0, 25:45, 40:60]).all()
      assert data.shape == labels.shape
      # Now test with more than two labels
      labels[55, 80] = 3
      full_prob_bf = random_walker(data, labels, beta=90, mode='bf',
                                   return_full_prob=True)
      assert (full_prob_bf[1, 25:45, 40:60] >=
              full_prob_bf[0, 25:45, 40:60]).all()
      assert len(full_prob_bf) == 3
      assert data.shape == labels.shape
  
  
  def test_2d_cg():
      lx = 70
      ly = 100
      data, labels = make_2d_syntheticdata(lx, ly)
      with expected_warnings(['"cg" mode' + '|' + SCIPY_EXPECTED]):
          labels_cg = random_walker(data, labels, beta=90, mode='cg')
      assert (labels_cg[25:45, 40:60] == 2).all()
      assert data.shape == labels.shape
      with expected_warnings(['"cg" mode' + '|' + SCIPY_EXPECTED]):
          full_prob = random_walker(data, labels, beta=90, mode='cg',
                                    return_full_prob=True)
      assert (full_prob[1, 25:45, 40:60] >=
              full_prob[0, 25:45, 40:60]).all()
      assert data.shape == labels.shape
      return data, labels_cg
  
  
  def test_2d_cg_mg():
      lx = 70
      ly = 100
      data, labels = make_2d_syntheticdata(lx, ly)
      expected = 'scipy.sparse.sparsetools|%s' % PYAMG_SCIPY_EXPECTED
      with expected_warnings([expected]):
          labels_cg_mg = random_walker(data, labels, beta=90, mode='cg_mg')
      assert (labels_cg_mg[25:45, 40:60] == 2).all()
      assert data.shape == labels.shape
      with expected_warnings([expected]):
          full_prob = random_walker(data, labels, beta=90, mode='cg_mg',
                                return_full_prob=True)
      assert (full_prob[1, 25:45, 40:60] >=
              full_prob[0, 25:45, 40:60]).all()
      assert data.shape == labels.shape
      return data, labels_cg_mg
  
  
  def test_types():
      lx = 70
      ly = 100
      data, labels = make_2d_syntheticdata(lx, ly)
      data = 255 * (data - data.min()) // (data.max() - data.min())
      data = data.astype(np.uint8)
      with expected_warnings([PYAMG_SCIPY_EXPECTED]):
          labels_cg_mg = random_walker(data, labels, beta=90, mode='cg_mg')
      assert (labels_cg_mg[25:45, 40:60] == 2).all()
      assert data.shape == labels.shape
      return data, labels_cg_mg
  
  
  def test_reorder_labels():
      lx = 70
      ly = 100
      data, labels = make_2d_syntheticdata(lx, ly)
      labels[labels == 2] = 4
      labels_bf = random_walker(data, labels, beta=90, mode='bf')
      assert (labels_bf[25:45, 40:60] == 2).all()
      assert data.shape == labels.shape
      return data, labels_bf
  
  
  def test_2d_inactive():
      lx = 70
      ly = 100
      data, labels = make_2d_syntheticdata(lx, ly)
      labels[10:20, 10:20] = -1
      labels[46:50, 33:38] = -2
      labels = random_walker(data, labels, beta=90)
      assert (labels.reshape((lx, ly))[25:45, 40:60] == 2).all()
      assert data.shape == labels.shape
      return data, labels
  
  
  def test_3d():
      n = 30
      lx, ly, lz = n, n, n
      data, labels = make_3d_syntheticdata(lx, ly, lz)
      with expected_warnings(['"cg" mode' + '|' + SCIPY_EXPECTED]):
          labels = random_walker(data, labels, mode='cg')
      assert (labels.reshape(data.shape)[13:17, 13:17, 13:17] == 2).all()
      assert data.shape == labels.shape
      return data, labels
  
  
  def test_3d_inactive():
      n = 30
      lx, ly, lz = n, n, n
      data, labels = make_3d_syntheticdata(lx, ly, lz)
      old_labels = np.copy(labels)
      labels[5:25, 26:29, 26:29] = -1
      after_labels = np.copy(labels)
      with expected_warnings(['"cg" mode|CObject type' + '|' + SCIPY_EXPECTED]):
          labels = random_walker(data, labels, mode='cg')
      assert (labels.reshape(data.shape)[13:17, 13:17, 13:17] == 2).all()
      assert data.shape == labels.shape
      return data, labels, old_labels, after_labels
  
  
  def test_multispectral_2d():
      lx, ly = 70, 100
      data, labels = make_2d_syntheticdata(lx, ly)
      data = data[..., np.newaxis].repeat(2, axis=-1)  # Expect identical output
      with expected_warnings(['"cg" mode' + '|' + SCIPY_EXPECTED]):
          multi_labels = random_walker(data, labels, mode='cg',
                                       multichannel=True)
      assert data[..., 0].shape == labels.shape
      with expected_warnings(['"cg" mode' + '|' + SCIPY_EXPECTED]):
          single_labels = random_walker(data[..., 0], labels, mode='cg')
      assert (multi_labels.reshape(labels.shape)[25:45, 40:60] == 2).all()
      assert data[..., 0].shape == labels.shape
      return data, multi_labels, single_labels, labels
  
  
  def test_multispectral_3d():
      n = 30
      lx, ly, lz = n, n, n
      data, labels = make_3d_syntheticdata(lx, ly, lz)
      data = data[..., np.newaxis].repeat(2, axis=-1)  # Expect identical output
      with expected_warnings(['"cg" mode' + '|' + SCIPY_EXPECTED]):
          multi_labels = random_walker(data, labels, mode='cg',
                                       multichannel=True)
      assert data[..., 0].shape == labels.shape
      with expected_warnings(['"cg" mode' + '|' + SCIPY_EXPECTED]):
          single_labels = random_walker(data[..., 0], labels, mode='cg')
      assert (multi_labels.reshape(labels.shape)[13:17, 13:17, 13:17] == 2).all()
      assert (single_labels.reshape(labels.shape)[13:17, 13:17, 13:17] == 2).all()
      assert data[..., 0].shape == labels.shape
      return data, multi_labels, single_labels, labels
  
  
  def test_spacing_0():
      n = 30
      lx, ly, lz = n, n, n
      data, _ = make_3d_syntheticdata(lx, ly, lz)
  
      # Rescale `data` along Z axis
      data_aniso = np.zeros((n, n, n // 2))
      for i, yz in enumerate(data):
          data_aniso[i, :, :] = resize(yz, (n, n // 2),
                                       mode='constant',
                                       anti_aliasing=False)
  
      # Generate new labels
      small_l = int(lx // 5)
      labels_aniso = np.zeros_like(data_aniso)
      labels_aniso[lx // 5, ly // 5, lz // 5] = 1
      labels_aniso[lx // 2 + small_l // 4,
                   ly // 2 - small_l // 4,
                   lz // 4 - small_l // 8] = 2
  
      # Test with `spacing` kwarg
      with expected_warnings(['"cg" mode' + '|' + SCIPY_EXPECTED]):
          labels_aniso = random_walker(data_aniso, labels_aniso, mode='cg',
                                       spacing=(1., 1., 0.5))
  
      assert (labels_aniso[13:17, 13:17, 7:9] == 2).all()
  
  
  def test_spacing_1():
      n = 30
      lx, ly, lz = n, n, n
      data, _ = make_3d_syntheticdata(lx, ly, lz)
  
      # Rescale `data` along Y axis
      # `resize` is not yet 3D capable, so this must be done by looping in 2D.
      data_aniso = np.zeros((n, n * 2, n))
      for i, yz in enumerate(data):
          data_aniso[i, :, :] = resize(yz, (n * 2, n),
                                       mode='constant',
                                       anti_aliasing=False)
  
      # Generate new labels
      small_l = int(lx // 5)
      labels_aniso = np.zeros_like(data_aniso)
      labels_aniso[lx // 5, ly // 5, lz // 5] = 1
      labels_aniso[lx // 2 + small_l // 4,
                   ly - small_l // 2,
                   lz // 2 - small_l // 4] = 2
  
      # Test with `spacing` kwarg
      # First, anisotropic along Y
      with expected_warnings(['"cg" mode' + '|' + SCIPY_EXPECTED]):
          labels_aniso = random_walker(data_aniso, labels_aniso, mode='cg',
                                       spacing=(1., 2., 1.))
      assert (labels_aniso[13:17, 26:34, 13:17] == 2).all()
  
      # Rescale `data` along X axis
      # `resize` is not yet 3D capable, so this must be done by looping in 2D.
      data_aniso = np.zeros((n, n * 2, n))
      for i in range(data.shape[1]):
          data_aniso[i, :, :] = resize(data[:, 1, :], (n * 2, n),
                                       mode='constant',
                                       anti_aliasing=False)
  
      # Generate new labels
      small_l = int(lx // 5)
      labels_aniso2 = np.zeros_like(data_aniso)
      labels_aniso2[lx // 5, ly // 5, lz // 5] = 1
      labels_aniso2[lx - small_l // 2,
                    ly // 2 + small_l // 4,
                    lz // 2 - small_l // 4] = 2
  
      # Anisotropic along X
      with expected_warnings(['"cg" mode' + '|' + SCIPY_EXPECTED]):
          labels_aniso2 = random_walker(data_aniso,
                                        labels_aniso2,
                                        mode='cg', spacing=(2., 1., 1.))
      assert (labels_aniso2[26:34, 13:17, 13:17] == 2).all()
  
  
  def test_trivial_cases():
      # When all voxels are labeled
      img = np.ones((10, 10))
      labels = np.ones((10, 10))
  
      with expected_warnings(["Returning provided labels"]):
          pass_through = random_walker(img, labels)
      np.testing.assert_array_equal(pass_through, labels)
  
      # When all voxels are labeled AND return_full_prob is True
      labels[:, :5] = 3
      expected = np.concatenate(((labels == 1)[..., np.newaxis],
                                 (labels == 3)[..., np.newaxis]), axis=2)
      with expected_warnings(["Returning provided labels"]):
          test = random_walker(img, labels, return_full_prob=True)
      np.testing.assert_array_equal(test, expected)
  
  
  def test_length2_spacing():
      # If this passes without raising an exception (warnings OK), the new
      #   spacing code is working properly.
      np.random.seed(42)
      img = np.ones((10, 10)) + 0.2 * np.random.normal(size=(10, 10))
      labels = np.zeros((10, 10), dtype=np.uint8)
      labels[2, 4] = 1
      labels[6, 8] = 4
      random_walker(img, labels, spacing=(1., 2.))
  
  
  def test_bad_inputs():
      # Too few dimensions
      img = np.ones(10)
      labels = np.arange(10)
      with testing.raises(ValueError):
          random_walker(img, labels)
      with testing.raises(ValueError):
          random_walker(img, labels, multichannel=True)
  
      # Too many dimensions
      np.random.seed(42)
      img = np.random.normal(size=(3, 3, 3, 3, 3))
      labels = np.arange(3 ** 5).reshape(img.shape)
      with testing.raises(ValueError):
          random_walker(img, labels)
      with testing.raises(ValueError):
          random_walker(img, labels, multichannel=True)
  
      # Spacing incorrect length
      img = np.random.normal(size=(10, 10))
      labels = np.zeros((10, 10))
      labels[2, 4] = 2
      labels[6, 8] = 5
      with testing.raises(ValueError):
          random_walker(img, labels, spacing=(1,))
  
      # Invalid mode
      img = np.random.normal(size=(10, 10))
      labels = np.zeros((10, 10))
      with testing.raises(ValueError):
          random_walker(img, labels, mode='bad')
  
  
  def test_isolated_seeds():
      np.random.seed(0)
      a = np.random.random((7, 7))
      mask = - np.ones(a.shape)
      # This pixel is an isolated seed
      mask[1, 1] = 1
      # Unlabeled pixels
      mask[3:, 3:] = 0
      # Seeds connected to unlabeled pixels
      mask[4, 4] = 2
      mask[6, 6] = 1
  
      # Test that no error is raised, and that labels of isolated seeds are OK
      res = random_walker(a, mask)
      assert res[1, 1] == 1
      res = random_walker(a, mask, return_full_prob=True)
      assert res[0, 1, 1] == 1
      assert res[1, 1, 1] == 0
   """Implementations restoration functions"""
  
  
  import numpy as np
  import numpy.random as npr
  from scipy.signal import fftconvolve, convolve
  
  from . import uft
  
  __keywords__ = "restoration, image, deconvolution"
  
  
  def wiener(image, psf, balance, reg=None, is_real=True, clip=True):
      """Wiener-Hunt deconvolution
  
      Return the deconvolution with a Wiener-Hunt approach (i.e. with
      Fourier diagonalisation).
  
      Parameters
      ----------
      image : (M, N) ndarray
         Input degraded image
      psf : ndarray
         Point Spread Function. This is assumed to be the impulse
         response (input image space) if the data-type is real, or the
         transfer function (Fourier space) if the data-type is
         complex. There is no constraints on the shape of the impulse
         response. The transfer function must be of shape `(M, N)` if
         `is_real is True`, `(M, N // 2 + 1)` otherwise (see
         `np.fft.rfftn`).
      balance : float
         The regularisation parameter value that tunes the balance
         between the data adequacy that improve frequency restoration
         and the prior adequacy that reduce frequency restoration (to
         avoid noise artifacts).
      reg : ndarray, optional
         The regularisation operator. The Laplacian by default. It can
         be an impulse response or a transfer function, as for the
         psf. Shape constraint is the same as for the `psf` parameter.
      is_real : boolean, optional
         True by default. Specify if ``psf`` and ``reg`` are provided
         with hermitian hypothesis, that is only half of the frequency
         plane is provided (due to the redundancy of Fourier transform
         of real signal). It's apply only if ``psf`` and/or ``reg`` are
         provided as transfer function.  For the hermitian property see
         ``uft`` module or ``np.fft.rfftn``.
      clip : boolean, optional
         True by default. If True, pixel values of the result above 1 or
         under -1 are thresholded for skimage pipeline compatibility.
  
      Returns
      -------
      im_deconv : (M, N) ndarray
         The deconvolved image.
  
      Examples
      --------
      >>> from skimage import color, data, restoration
      >>> img = color.rgb2gray(data.astronaut())
      >>> from scipy.signal import convolve2d
      >>> psf = np.ones((5, 5)) / 25
      >>> img = convolve2d(img, psf, 'same')
      >>> img += 0.1 * img.std() * np.random.standard_normal(img.shape)
      >>> deconvolved_img = restoration.wiener(img, psf, 1100)
  
      Notes
      -----
      This function applies the Wiener filter to a noisy and degraded
      image by an impulse response (or PSF). If the data model is
  
      .. math:: y = Hx + n
  
      where :math:`n` is noise, :math:`H` the PSF and :math:`x` the
      unknown original image, the Wiener filter is
  
      .. math::
         \hat x = F^\dagger (|\Lambda_H|^2 + \lambda |\Lambda_D|^2)
         \Lambda_H^\dagger F y
  
      where :math:`F` and :math:`F^\dagger` are the Fourier and inverse
      Fourier transfroms respectively, :math:`\Lambda_H` the transfer
      function (or the Fourier transfrom of the PSF, see [Hunt] below)
      and :math:`\Lambda_D` the filter to penalize the restored image
      frequencies (Laplacian by default, that is penalization of high
      frequency). The parameter :math:`\lambda` tunes the balance
      between the data (that tends to increase high frequency, even
      those coming from noise), and the regularization.
  
      These methods are then specific to a prior model. Consequently,
      the application or the true image nature must corresponds to the
      prior model. By default, the prior model (Laplacian) introduce
      image smoothness or pixel correlation. It can also be interpreted
      as high-frequency penalization to compensate the instability of
      the solution with respect to the data (sometimes called noise
      amplification or "explosive" solution).
  
      Finally, the use of Fourier space implies a circulant property of
      :math:`H`, see [Hunt].
  
      References
      ----------
      .. [1] François Orieux, Jean-François Giovannelli, and Thomas
             Rodet, "Bayesian estimation of regularization and point
             spread function parameters for Wiener-Hunt deconvolution",
             J. Opt. Soc. Am. A 27, 1593-1607 (2010)
  
             http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-27-7-1593
  
             http://research.orieux.fr/files/papers/OGR-JOSA10.pdf
  
      .. [2] B. R. Hunt "A matrix theory proof of the discrete
             convolution theorem", IEEE Trans. on Audio and
             Electroacoustics, vol. au-19, no. 4, pp. 285-288, dec. 1971
      """
      if reg is None:
          reg, _ = uft.laplacian(image.ndim, image.shape, is_real=is_real)
      if not np.iscomplexobj(reg):
          reg = uft.ir2tf(reg, image.shape, is_real=is_real)
  
      if psf.shape != reg.shape:
          trans_func = uft.ir2tf(psf, image.shape, is_real=is_real)
      else:
          trans_func = psf
  
      wiener_filter = np.conj(trans_func) / (np.abs(trans_func) ** 2 +
                                             balance * np.abs(reg) ** 2)
      if is_real:
          deconv = uft.uirfft2(wiener_filter * uft.urfft2(image),
                               shape=image.shape)
      else:
          deconv = uft.uifft2(wiener_filter * uft.ufft2(image))
  
      if clip:
          deconv[deconv > 1] = 1
          deconv[deconv < -1] = -1
  
      return deconv
  
  
  def unsupervised_wiener(image, psf, reg=None, user_params=None, is_real=True,
                          clip=True):
      """Unsupervised Wiener-Hunt deconvolution.
  
      Return the deconvolution with a Wiener-Hunt approach, where the
      hyperparameters are automatically estimated. The algorithm is a
      stochastic iterative process (Gibbs sampler) described in the
      reference below. See also ``wiener`` function.
  
      Parameters
      ----------
      image : (M, N) ndarray
         The input degraded image.
      psf : ndarray
         The impulse response (input image's space) or the transfer
         function (Fourier space). Both are accepted. The transfer
         function is automatically recognized as being complex
         (``np.iscomplexobj(psf)``).
      reg : ndarray, optional
         The regularisation operator. The Laplacian by default. It can
         be an impulse response or a transfer function, as for the psf.
      user_params : dict
         Dictionary of parameters for the Gibbs sampler. See below.
      clip : boolean, optional
         True by default. If true, pixel values of the result above 1 or
         under -1 are thresholded for skimage pipeline compatibility.
  
      Returns
      -------
      x_postmean : (M, N) ndarray
         The deconvolved image (the posterior mean).
      chains : dict
         The keys ``noise`` and ``prior`` contain the chain list of
         noise and prior precision respectively.
  
      Other parameters
      ----------------
      The keys of ``user_params`` are:
  
      threshold : float
         The stopping criterion: the norm of the difference between to
         successive approximated solution (empirical mean of object
         samples, see Notes section). 1e-4 by default.
      burnin : int
         The number of sample to ignore to start computation of the
         mean. 15 by default.
      min_iter : int
         The minimum number of iterations. 30 by default.
      max_iter : int
         The maximum number of iterations if ``threshold`` is not
         satisfied. 200 by default.
      callback : callable (None by default)
         A user provided callable to which is passed, if the function
         exists, the current image sample for whatever purpose. The user
         can store the sample, or compute other moments than the
         mean. It has no influence on the algorithm execution and is
         only for inspection.
  
      Examples
      --------
      >>> from skimage import color, data, restoration
      >>> img = color.rgb2gray(data.astronaut())
      >>> from scipy.signal import convolve2d
      >>> psf = np.ones((5, 5)) / 25
      >>> img = convolve2d(img, psf, 'same')
      >>> img += 0.1 * img.std() * np.random.standard_normal(img.shape)
      >>> deconvolved_img = restoration.unsupervised_wiener(img, psf)
  
      Notes
      -----
      The estimated image is design as the posterior mean of a
      probability law (from a Bayesian analysis). The mean is defined as
      a sum over all the possible images weighted by their respective
      probability. Given the size of the problem, the exact sum is not
      tractable. This algorithm use of MCMC to draw image under the
      posterior law. The practical idea is to only draw highly probable
      images since they have the biggest contribution to the mean. At the
      opposite, the less probable images are drawn less often since
      their contribution is low. Finally the empirical mean of these
      samples give us an estimation of the mean, and an exact
      computation with an infinite sample set.
  
      References
      ----------
      .. [1] François Orieux, Jean-François Giovannelli, and Thomas
             Rodet, "Bayesian estimation of regularization and point
             spread function parameters for Wiener-Hunt deconvolution",
             J. Opt. Soc. Am. A 27, 1593-1607 (2010)
  
             http://www.opticsinfobase.org/josaa/abstract.cfm?URI=josaa-27-7-1593
  
             http://research.orieux.fr/files/papers/OGR-JOSA10.pdf
      """
      params = {'threshold': 1e-4, 'max_iter': 200,
                'min_iter': 30, 'burnin': 15, 'callback': None}
      params.update(user_params or {})
  
      if reg is None:
          reg, _ = uft.laplacian(image.ndim, image.shape, is_real=is_real)
      if not np.iscomplexobj(reg):
          reg = uft.ir2tf(reg, image.shape, is_real=is_real)
  
      if psf.shape != reg.shape:
          trans_fct = uft.ir2tf(psf, image.shape,  is_real=is_real)
      else:
          trans_fct = psf
  
      # The mean of the object
      x_postmean = np.zeros(trans_fct.shape)
      # The previous computed mean in the iterative loop
      prev_x_postmean = np.zeros(trans_fct.shape)
  
      # Difference between two successive mean
      delta = np.NAN
  
      # Initial state of the chain
      gn_chain, gx_chain = [1], [1]
  
      # The correlation of the object in Fourier space (if size is big,
      # this can reduce computation time in the loop)
      areg2 = np.abs(reg) ** 2
      atf2 = np.abs(trans_fct) ** 2
  
      # The Fourier transfrom may change the image.size attribut, so we
      # store it.
      if is_real:
          data_spectrum = uft.urfft2(image.astype(np.float))
      else:
          data_spectrum = uft.ufft2(image.astype(np.float))
  
      # Gibbs sampling
      for iteration in range(params['max_iter']):
          # Sample of Eq. 27 p(circX^k | gn^k-1, gx^k-1, y).
  
          # weighting (correlation in direct space)
          precision = gn_chain[-1] * atf2 + gx_chain[-1] * areg2  # Eq. 29
          excursion = np.sqrt(0.5) / np.sqrt(precision) * (
              np.random.standard_normal(data_spectrum.shape) +
              1j * np.random.standard_normal(data_spectrum.shape))
  
          # mean Eq. 30 (RLS for fixed gn, gamma0 and gamma1 ...)
          wiener_filter = gn_chain[-1] * np.conj(trans_fct) / precision
  
          # sample of X in Fourier space
          x_sample = wiener_filter * data_spectrum + excursion
          if params['callback']:
              params['callback'](x_sample)
  
          # sample of Eq. 31 p(gn | x^k, gx^k, y)
          gn_chain.append(npr.gamma(image.size / 2,
                                    2 / uft.image_quad_norm(data_spectrum -
                                                            x_sample *
                                                            trans_fct)))
  
          # sample of Eq. 31 p(gx | x^k, gn^k-1, y)
          gx_chain.append(npr.gamma((image.size - 1) / 2,
                                    2 / uft.image_quad_norm(x_sample * reg)))
  
          # current empirical average
          if iteration > params['burnin']:
              x_postmean = prev_x_postmean + x_sample
  
          if iteration > (params['burnin'] + 1):
              current = x_postmean / (iteration - params['burnin'])
              previous = prev_x_postmean / (iteration - params['burnin'] - 1)
  
              delta = np.sum(np.abs(current - previous)) / \
                  np.sum(np.abs(x_postmean)) / (iteration - params['burnin'])
  
          prev_x_postmean = x_postmean
  
          # stop of the algorithm
          if (iteration > params['min_iter']) and (delta < params['threshold']):
              break
  
      # Empirical average \approx POSTMEAN Eq. 44
      x_postmean = x_postmean / (iteration - params['burnin'])
      if is_real:
          x_postmean = uft.uirfft2(x_postmean, shape=image.shape)
      else:
          x_postmean = uft.uifft2(x_postmean)
  
      if clip:
          x_postmean[x_postmean > 1] = 1
          x_postmean[x_postmean < -1] = -1
  
      return (x_postmean, {'noise': gn_chain, 'prior': gx_chain})
  
  
  def richardson_lucy(image, psf, iterations=50, clip=True):
      """Richardson-Lucy deconvolution.
  
      Parameters
      ----------
      image : ndarray
         Input degraded image (can be N dimensional).
      psf : ndarray
         The point spread function.
      iterations : int
         Number of iterations. This parameter plays the role of
         regularisation.
      clip : boolean, optional
         True by default. If true, pixel value of the result above 1 or
         under -1 are thresholded for skimage pipeline compatibility.
  
      Returns
      -------
      im_deconv : ndarray
         The deconvolved image.
  
      Examples
      --------
      >>> from skimage import color, data, restoration
      >>> camera = color.rgb2gray(data.camera())
      >>> from scipy.signal import convolve2d
      >>> psf = np.ones((5, 5)) / 25
      >>> camera = convolve2d(camera, psf, 'same')
      >>> camera += 0.1 * camera.std() * np.random.standard_normal(camera.shape)
      >>> deconvolved = restoration.richardson_lucy(camera, psf, 5)
  
      References
      ----------
      .. [1] http://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution
      """
      # compute the times for direct convolution and the fft method. The fft is of
      # complexity O(N log(N)) for each dimension and the direct method does
      # straight arithmetic (and is O(n*k) to add n elements k times)
      direct_time = np.prod(image.shape + psf.shape)
      fft_time =  np.sum([n*np.log(n) for n in image.shape + psf.shape])
  
      # see whether the fourier transform convolution method or the direct
      # convolution method is faster (discussed in scikit-image PR #1792)
      time_ratio = 40.032 * fft_time / direct_time
  
      if time_ratio <= 1 or len(image.shape) > 2:
          convolve_method = fftconvolve
      else:
          convolve_method = convolve
  
      image = image.astype(np.float)
      psf = psf.astype(np.float)
      im_deconv = 0.5 * np.ones(image.shape)
      psf_mirror = psf[::-1, ::-1]
  
      for _ in range(iterations):
          relative_blur = image / convolve_method(im_deconv, psf, 'same')
          im_deconv *= convolve_method(relative_blur, psf_mirror, 'same')
  
      if clip:
          im_deconv[im_deconv > 1] = 1
          im_deconv[im_deconv < -1] = -1
  
      return im_deconv
   import numpy as np
  import scipy.ndimage as ndi
  from ..segmentation import relabel_sequential
  from .. import measure
  from ..filters import rank_order
  
  
  def _get_high_intensity_peaks(image, mask, num_peaks):
      """
      Return the highest intensity peak coordinates.
      """
      # get coordinates of peaks
      coord = np.nonzero(mask)
      # select num_peaks peaks
      if len(coord[0]) > num_peaks:
          intensities = image[coord]
          idx_maxsort = np.argsort(intensities)
          coord = np.transpose(coord)[idx_maxsort][-num_peaks:]
      else:
          coord = np.column_stack(coord)
      # Higest peak first
      return coord[::-1]
  
  
  def peak_local_max(image, min_distance=1, threshold_abs=None,
                     threshold_rel=None, exclude_border=True, indices=True,
                     num_peaks=np.inf, footprint=None, labels=None,
                     num_peaks_per_label=np.inf):
      """Find peaks in an image as coordinate list or boolean mask.
  
      Peaks are the local maxima in a region of `2 * min_distance + 1`
      (i.e. peaks are separated by at least `min_distance`).
  
      If peaks are flat (i.e. multiple adjacent pixels have identical
      intensities), the coordinates of all such pixels are returned.
  
      If both `threshold_abs` and `threshold_rel` are provided, the maximum
      of the two is chosen as the minimum intensity threshold of peaks.
  
      Parameters
      ----------
      image : ndarray
          Input image.
      min_distance : int, optional
          Minimum number of pixels separating peaks in a region of `2 *
          min_distance + 1` (i.e. peaks are separated by at least
          `min_distance`).
          To find the maximum number of peaks, use `min_distance=1`.
      threshold_abs : float, optional
          Minimum intensity of peaks. By default, the absolute threshold is
          the minimum intensity of the image.
      threshold_rel : float, optional
          Minimum intensity of peaks, calculated as `max(image) * threshold_rel`.
      exclude_border : int, optional
          If nonzero, `exclude_border` excludes peaks from
          within `exclude_border`-pixels of the border of the image.
      indices : bool, optional
          If True, the output will be an array representing peak
          coordinates.  If False, the output will be a boolean array shaped as
          `image.shape` with peaks present at True elements.
      num_peaks : int, optional
          Maximum number of peaks. When the number of peaks exceeds `num_peaks`,
          return `num_peaks` peaks based on highest peak intensity.
      footprint : ndarray of bools, optional
          If provided, `footprint == 1` represents the local region within which
          to search for peaks at every point in `image`.  Overrides
          `min_distance` (also for `exclude_border`).
      labels : ndarray of ints, optional
          If provided, each unique region `labels == value` represents a unique
          region to search for peaks. Zero is reserved for background.
      num_peaks_per_label : int, optional
          Maximum number of peaks for each label.
  
      Returns
      -------
      output : ndarray or ndarray of bools
  
          * If `indices = True`  : (row, column, ...) coordinates of peaks.
          * If `indices = False` : Boolean array shaped like `image`, with peaks
            represented by True values.
  
      Notes
      -----
      The peak local maximum function returns the coordinates of local peaks
      (maxima) in an image. A maximum filter is used for finding local maxima.
      This operation dilates the original image. After comparison of the dilated
      and original image, this function returns the coordinates or a mask of the
      peaks where the dilated image equals the original image.
  
      Examples
      --------
      >>> img1 = np.zeros((7, 7))
      >>> img1[3, 4] = 1
      >>> img1[3, 2] = 1.5
      >>> img1
      array([[ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
             [ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
             [ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
             [ 0. ,  0. ,  1.5,  0. ,  1. ,  0. ,  0. ],
             [ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
             [ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
             [ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ]])
  
      >>> peak_local_max(img1, min_distance=1)
      array([[3, 4],
             [3, 2]])
  
      >>> peak_local_max(img1, min_distance=2)
      array([[3, 2]])
  
      >>> img2 = np.zeros((20, 20, 20))
      >>> img2[10, 10, 10] = 1
      >>> peak_local_max(img2, exclude_border=0)
      array([[10, 10, 10]])
  
      """
      if type(exclude_border) == bool:
          exclude_border = min_distance if exclude_border else 0
  
      out = np.zeros_like(image, dtype=np.bool)
  
      # In the case of labels, recursively build and return an output
      # operating on each label separately
      if labels is not None:
          label_values = np.unique(labels)
          # Reorder label values to have consecutive integers (no gaps)
          if np.any(np.diff(label_values) != 1):
              mask = labels >= 1
              labels[mask] = 1 + rank_order(labels[mask])[0].astype(labels.dtype)
          labels = labels.astype(np.int32)
  
          # New values for new ordering
          label_values = np.unique(labels)
          for label in label_values[label_values != 0]:
              maskim = (labels == label)
              out += peak_local_max(image * maskim, min_distance=min_distance,
                                    threshold_abs=threshold_abs,
                                    threshold_rel=threshold_rel,
                                    exclude_border=exclude_border,
                                    indices=False, num_peaks=num_peaks_per_label,
                                    footprint=footprint, labels=None)
  
          # Select highest intensities (num_peaks)
          coordinates = _get_high_intensity_peaks(image, out, num_peaks)
  
          if indices is True:
              return coordinates
          else:
              nd_indices = tuple(coordinates.T)
              out[nd_indices] = True
              return out
  
      if np.all(image == image.flat[0]):
          if indices is True:
              return np.empty((0, 2), np.int)
          else:
              return out
  
      # Non maximum filter
      if footprint is not None:
          image_max = ndi.maximum_filter(image, footprint=footprint,
                                         mode='constant')
      else:
          size = 2 * min_distance + 1
          image_max = ndi.maximum_filter(image, size=size, mode='constant')
      mask = image == image_max
  
      if exclude_border:
          # zero out the image borders
          for i in range(mask.ndim):
              mask = mask.swapaxes(0, i)
              remove = (footprint.shape[i] if footprint is not None
                        else 2 * exclude_border)
              mask[:remove // 2] = mask[-remove // 2:] = False
              mask = mask.swapaxes(0, i)
  
      # find top peak candidates above a threshold
      thresholds = []
      if threshold_abs is None:
          threshold_abs = image.min()
      thresholds.append(threshold_abs)
      if threshold_rel is not None:
          thresholds.append(threshold_rel * image.max())
      if thresholds:
          mask &= image > max(thresholds)
  
      # Select highest intensities (num_peaks)
      coordinates = _get_high_intensity_peaks(image, mask, num_peaks)
  
      if indices is True:
          return coordinates
      else:
          nd_indices = tuple(coordinates.T)
          out[nd_indices] = True
          return out
  
  
  def _prominent_peaks(image, min_xdistance=1, min_ydistance=1,
                       threshold=None, num_peaks=np.inf):
      """Return peaks with non-maximum suppression.
  
      Identifies most prominent features separated by certain distances.
      Non-maximum suppression with different sizes is applied separately
      in the first and second dimension of the image to identify peaks.
  
      Parameters
      ----------
      image : (M, N) ndarray
          Input image.
      min_xdistance : int
          Minimum distance separating features in the x dimension.
      min_ydistance : int
          Minimum distance separating features in the y dimension.
      threshold : float
          Minimum intensity of peaks. Default is `0.5 * max(image)`.
      num_peaks : int
          Maximum number of peaks. When the number of peaks exceeds `num_peaks`,
          return `num_peaks` coordinates based on peak intensity.
  
      Returns
      -------
      intensity, xcoords, ycoords : tuple of array
          Peak intensity values, x and y indices.
      """
  
      img = image.copy()
      rows, cols = img.shape
  
      if threshold is None:
          threshold = 0.5 * np.max(img)
  
      ycoords_size = 2 * min_ydistance + 1
      xcoords_size = 2 * min_xdistance + 1
      img_max = ndi.maximum_filter1d(img, size=ycoords_size, axis=0,
                                     mode='constant', cval=0)
      img_max = ndi.maximum_filter1d(img_max, size=xcoords_size, axis=1,
                                     mode='constant', cval=0)
      mask = (img == img_max)
      img *= mask
      img_t = img > threshold
  
      label_img = measure.label(img_t)
      props = measure.regionprops(label_img, img_max)
  
      # Sort the list of peaks by intensity, not left-right, so larger peaks
      # in Hough space cannot be arbitrarily suppressed by smaller neighbors
      props = sorted(props, key=lambda x: x.max_intensity)[::-1]
      coords = np.array([np.round(p.centroid) for p in props], dtype=int)
  
      img_peaks = []
      ycoords_peaks = []
      xcoords_peaks = []
  
      # relative coordinate grid for local neighbourhood suppression
      ycoords_ext, xcoords_ext = np.mgrid[-min_ydistance:min_ydistance + 1,
                                          -min_xdistance:min_xdistance + 1]
  
      for ycoords_idx, xcoords_idx in coords:
          accum = img_max[ycoords_idx, xcoords_idx]
          if accum > threshold:
              # absolute coordinate grid for local neighbourhood suppression
              ycoords_nh = ycoords_idx + ycoords_ext
              xcoords_nh = xcoords_idx + xcoords_ext
  
              # no reflection for distance neighbourhood
              ycoords_in = np.logical_and(ycoords_nh > 0, ycoords_nh < rows)
              ycoords_nh = ycoords_nh[ycoords_in]
              xcoords_nh = xcoords_nh[ycoords_in]
  
              # reflect xcoords and assume xcoords are continuous,
              # e.g. for angles:
              # (..., 88, 89, -90, -89, ..., 89, -90, -89, ...)
              xcoords_low = xcoords_nh < 0
              ycoords_nh[xcoords_low] = rows - ycoords_nh[xcoords_low]
              xcoords_nh[xcoords_low] += cols
              xcoords_high = xcoords_nh >= cols
              ycoords_nh[xcoords_high] = rows - ycoords_nh[xcoords_high]
              xcoords_nh[xcoords_high] -= cols
  
              # suppress neighbourhood
              img_max[ycoords_nh, xcoords_nh] = 0
  
              # add current feature to peaks
              img_peaks.append(accum)
              ycoords_peaks.append(ycoords_idx)
              xcoords_peaks.append(xcoords_idx)
  
      img_peaks = np.array(img_peaks)
      ycoords_peaks = np.array(ycoords_peaks)
      xcoords_peaks = np.array(xcoords_peaks)
  
      if num_peaks < len(img_peaks):
          idx_maxsort = np.argsort(img_peaks)[::-1][:num_peaks]
          img_peaks = img_peaks[idx_maxsort]
          ycoords_peaks = ycoords_peaks[idx_maxsort]
          xcoords_peaks = xcoords_peaks[idx_maxsort]
  
      return img_peaks, xcoords_peaks, ycoords_peaks
  
apis: numpy.zeros_like numpy.ones numpy.zeros numpy.random.seed numpy.testing.assert_array_equal numpy.copy numpy.random.randn numpy.arange numpy.random.random numpy.max numpy.random.normal numpy.concatenate
apis: numpy.ones numpy.random.standard_normal numpy.zeros numpy.conj numpy.abs numpy.log numpy.prod numpy.sqrt numpy.iscomplexobj
apis: numpy.zeros_like numpy.transpose numpy.empty scipy.ndimage.maximum_filter numpy.round numpy.diff numpy.logical_and numpy.argsort numpy.column_stack numpy.all numpy.max scipy.ndimage.maximum_filter1d numpy.array numpy.nonzero numpy.unique