code: # -*- coding: utf-8 -*-
  """
  radon.py - Radon and inverse radon transforms
  
  Based on code of Justin K. Romberg
  (http://www.clear.rice.edu/elec431/projects96/DSP/bpanalysis.html)
  J. Gillam and Chris Griffin.
  
  References:
      -B.R. Ramesh, N. Srinivasa, K. Rajgopal, "An Algorithm for Computing
      the Discrete Radon Transform With Some Applications", Proceedings of
      the Fourth IEEE Region 10 International Conference, TENCON '89, 1989.
      -A. C. Kak, Malcolm Slaney, "Principles of Computerized Tomographic
      Imaging", IEEE Press 1988.
  """
  from __future__ import division
  import numpy as np
  from scipy.fftpack import fft, ifft, fftfreq
  from scipy.interpolate import interp1d
  from ._warps_cy import _warp_fast
  from ._radon_transform import sart_projection_update
  from .. import util
  from warnings import warn
  
  
  __all__ = ["radon", "iradon", "iradon_sart"]
  
  
  def radon(image, theta=None, circle=False):
      """
      Calculates the radon transform of an image given specified
      projection angles.
  
      Parameters
      ----------
      image : array_like, dtype=float
          Input image. The rotation axis will be located in the pixel with
          indices ``(image.shape[0] // 2, image.shape[1] // 2)``.
      theta : array_like, dtype=float, optional (default np.arange(180))
          Projection angles (in degrees).
      circle : boolean, optional
          Assume image is zero outside the inscribed circle, making the
          width of each projection (the first dimension of the sinogram)
          equal to ``min(image.shape)``.
  
      Returns
      -------
      radon_image : ndarray
          Radon transform (sinogram).  The tomography rotation axis will lie
          at the pixel index ``radon_image.shape[0] // 2`` along the 0th
          dimension of ``radon_image``.
  
      """
      if image.ndim != 2:
          raise ValueError('The input image must be 2-D')
      if theta is None:
          theta = np.arange(180)
  
      if circle:
          radius = min(image.shape) // 2
          c0, c1 = np.ogrid[0:image.shape[0], 0:image.shape[1]]
          reconstruction_circle = ((c0 - image.shape[0] // 2) ** 2
                                   + (c1 - image.shape[1] // 2) ** 2)
          reconstruction_circle = reconstruction_circle <= radius ** 2
          if not np.all(reconstruction_circle | (image == 0)):
              warn('Radon transform: image must be zero outside the '
                   'reconstruction circle')
          # Crop image to make it square
          slices = []
          for d in (0, 1):
              if image.shape[d] > min(image.shape):
                  excess = image.shape[d] - min(image.shape)
                  slices.append(slice(int(np.ceil(excess / 2)),
                                      int(np.ceil(excess / 2)
                                          + min(image.shape))))
              else:
                  slices.append(slice(None))
          slices = tuple(slices)
          padded_image = image[slices]
      else:
          diagonal = np.sqrt(2) * max(image.shape)
          pad = [int(np.ceil(diagonal - s)) for s in image.shape]
          new_center = [(s + p) // 2 for s, p in zip(image.shape, pad)]
          old_center = [s // 2 for s in image.shape]
          pad_before = [nc - oc for oc, nc in zip(old_center, new_center)]
          pad_width = [(pb, p - pb) for pb, p in zip(pad_before, pad)]
          padded_image = util.pad(image, pad_width, mode='constant',
                                  constant_values=0)
      # padded_image is always square
      assert padded_image.shape[0] == padded_image.shape[1]
      radon_image = np.zeros((padded_image.shape[0], len(theta)))
      center = padded_image.shape[0] // 2
  
      shift0 = np.array([[1, 0, -center],
                         [0, 1, -center],
                         [0, 0, 1]])
      shift1 = np.array([[1, 0, center],
                         [0, 1, center],
                         [0, 0, 1]])
  
      def build_rotation(theta):
          T = np.deg2rad(theta)
          R = np.array([[np.cos(T), np.sin(T), 0],
                        [-np.sin(T), np.cos(T), 0],
                        [0, 0, 1]])
          return shift1.dot(R).dot(shift0)
  
      for i in range(len(theta)):
          rotated = _warp_fast(padded_image, build_rotation(theta[i]))
          radon_image[:, i] = rotated.sum(0)
      return radon_image
  
  
  def _sinogram_circle_to_square(sinogram):
      diagonal = int(np.ceil(np.sqrt(2) * sinogram.shape[0]))
      pad = diagonal - sinogram.shape[0]
      old_center = sinogram.shape[0] // 2
      new_center = diagonal // 2
      pad_before = new_center - old_center
      pad_width = ((pad_before, pad - pad_before), (0, 0))
      return util.pad(sinogram, pad_width, mode='constant', constant_values=0)
  
  
  def iradon(radon_image, theta=None, output_size=None,
             filter="ramp", interpolation="linear", circle=False):
      """
      Inverse radon transform.
  
      Reconstruct an image from the radon transform, using the filtered
      back projection algorithm.
  
      Parameters
      ----------
      radon_image : array_like, dtype=float
          Image containing radon transform (sinogram). Each column of
          the image corresponds to a projection along a different angle. The
          tomography rotation axis should lie at the pixel index
          ``radon_image.shape[0] // 2`` along the 0th dimension of
          ``radon_image``.
      theta : array_like, dtype=float, optional
          Reconstruction angles (in degrees). Default: m angles evenly spaced
          between 0 and 180 (if the shape of `radon_image` is (N, M)).
      output_size : int
          Number of rows and columns in the reconstruction.
      filter : str, optional (default ramp)
          Filter used in frequency domain filtering. Ramp filter used by default.
          Filters available: ramp, shepp-logan, cosine, hamming, hann.
          Assign None to use no filter.
      interpolation : str, optional (default 'linear')
          Interpolation method used in reconstruction. Methods available:
          'linear', 'nearest', and 'cubic' ('cubic' is slow).
      circle : boolean, optional
          Assume the reconstructed image is zero outside the inscribed circle.
          Also changes the default output_size to match the behaviour of
          ``radon`` called with ``circle=True``.
  
      Returns
      -------
      reconstructed : ndarray
          Reconstructed image. The rotation axis will be located in the pixel
          with indices
          ``(reconstructed.shape[0] // 2, reconstructed.shape[1] // 2)``.
  
      Notes
      -----
      It applies the Fourier slice theorem to reconstruct an image by
      multiplying the frequency domain of the filter with the FFT of the
      projection data. This algorithm is called filtered back projection.
  
      """
      if radon_image.ndim != 2:
          raise ValueError('The input image must be 2-D')
      if theta is None:
          m, n = radon_image.shape
          theta = np.linspace(0, 180, n, endpoint=False)
      else:
          theta = np.asarray(theta)
      if len(theta) != radon_image.shape[1]:
          raise ValueError("The given ``theta`` does not match the number of "
                           "projections in ``radon_image``.")
      interpolation_types = ('linear', 'nearest', 'cubic')
      if not interpolation in interpolation_types:
          raise ValueError("Unknown interpolation: %s" % interpolation)
      if not output_size:
          # If output size not specified, estimate from input radon image
          if circle:
              output_size = radon_image.shape[0]
          else:
              output_size = int(np.floor(np.sqrt((radon_image.shape[0]) ** 2
                                                 / 2.0)))
      if circle:
          radon_image = _sinogram_circle_to_square(radon_image)
  
      th = (np.pi / 180.0) * theta
      # resize image to next power of two (but no less than 64) for
      # Fourier analysis; speeds up Fourier and lessens artifacts
      projection_size_padded = \
          max(64, int(2 ** np.ceil(np.log2(2 * radon_image.shape[0]))))
      pad_width = ((0, projection_size_padded - radon_image.shape[0]), (0, 0))
      img = util.pad(radon_image, pad_width, mode='constant', constant_values=0)
  
      # Construct the Fourier filter
      f = fftfreq(projection_size_padded).reshape(-1, 1)   # digital frequency
      omega = 2 * np.pi * f                                # angular frequency
      fourier_filter = 2 * np.abs(f)                       # ramp filter
      if filter == "ramp":
          pass
      elif filter == "shepp-logan":
          # Start from first element to avoid divide by zero
          fourier_filter[1:] = fourier_filter[1:] * np.sin(omega[1:]) / omega[1:]
      elif filter == "cosine":
          fourier_filter *= np.cos(omega)
      elif filter == "hamming":
          fourier_filter *= (0.54 + 0.46 * np.cos(omega / 2))
      elif filter == "hann":
          fourier_filter *= (1 + np.cos(omega / 2)) / 2
      elif filter is None:
          fourier_filter[:] = 1
      else:
          raise ValueError("Unknown filter: %s" % filter)
      # Apply filter in Fourier domain
      projection = fft(img, axis=0) * fourier_filter
      radon_filtered = np.real(ifft(projection, axis=0))
  
      # Resize filtered image back to original size
      radon_filtered = radon_filtered[:radon_image.shape[0], :]
      reconstructed = np.zeros((output_size, output_size))
      # Determine the center of the projections (= center of sinogram)
      mid_index = radon_image.shape[0] // 2
  
      [X, Y] = np.mgrid[0:output_size, 0:output_size]
      xpr = X - int(output_size) // 2
      ypr = Y - int(output_size) // 2
  
      # Reconstruct image by interpolation
      for i in range(len(theta)):
          t = ypr * np.cos(th[i]) - xpr * np.sin(th[i])
          x = np.arange(radon_filtered.shape[0]) - mid_index
          if interpolation == 'linear':
              backprojected = np.interp(t, x, radon_filtered[:, i],
                                        left=0, right=0)
          else:
              interpolant = interp1d(x, radon_filtered[:, i], kind=interpolation,
                                     bounds_error=False, fill_value=0)
              backprojected = interpolant(t)
          reconstructed += backprojected
      if circle:
          radius = output_size // 2
          reconstruction_circle = (xpr ** 2 + ypr ** 2) <= radius ** 2
          reconstructed[~reconstruction_circle] = 0.
  
      return reconstructed * np.pi / (2 * len(th))
  
  
  def order_angles_golden_ratio(theta):
      """
      Order angles to reduce the amount of correlated information
      in subsequent projections.
  
      Parameters
      ----------
      theta : 1D array of floats
          Projection angles in degrees. Duplicate angles are not allowed.
  
      Returns
      -------
      indices_generator : generator yielding unsigned integers
          The returned generator yields indices into ``theta`` such that
          ``theta[indices]`` gives the approximate golden ratio ordering
          of the projections. In total, ``len(theta)`` indices are yielded.
          All non-negative integers < ``len(theta)`` are yielded exactly once.
  
      Notes
      -----
      The method used here is that of the golden ratio introduced
      by T. Kohler.
  
      References
      ----------
      .. [1] Kohler, T. "A projection access scheme for iterative
             reconstruction based on the golden section." Nuclear Science
             Symposium Conference Record, 2004 IEEE. Vol. 6. IEEE, 2004.
      .. [2] Winkelmann, Stefanie, et al. "An optimal radial profile order
             based on the Golden Ratio for time-resolved MRI."
             Medical Imaging, IEEE Transactions on 26.1 (2007): 68-76.
      """
      interval = 180
  
      def angle_distance(a, b):
          difference = a - b
          return min(abs(difference % interval), abs(difference % -interval))
  
      remaining = list(np.argsort(theta))   # indices into theta
      # yield an arbitrary angle to start things off
      index = remaining.pop(0)
      angle = theta[index]
      yield index
      # determine subsequent angles using the golden ratio method
      angle_increment = interval * (1 - (np.sqrt(5) - 1) / 2)
      while remaining:
          angle = (angle + angle_increment) % interval
          insert_point = np.searchsorted(theta[remaining], angle)
          index_below = insert_point - 1
          index_above = 0 if insert_point == len(remaining) else insert_point
          distance_below = angle_distance(angle, theta[remaining[index_below]])
          distance_above = angle_distance(angle, theta[remaining[index_above]])
          if distance_below < distance_above:
              yield remaining.pop(index_below)
          else:
              yield remaining.pop(index_above)
  
  
  def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
                  clip=None, relaxation=0.15):
      """
      Inverse radon transform
  
      Reconstruct an image from the radon transform, using a single iteration of
      the Simultaneous Algebraic Reconstruction Technique (SART) algorithm.
  
      Parameters
      ----------
      radon_image : 2D array, dtype=float
          Image containing radon transform (sinogram). Each column of
          the image corresponds to a projection along a different angle. The
          tomography rotation axis should lie at the pixel index
          ``radon_image.shape[0] // 2`` along the 0th dimension of
          ``radon_image``.
      theta : 1D array, dtype=float, optional
          Reconstruction angles (in degrees). Default: m angles evenly spaced
          between 0 and 180 (if the shape of `radon_image` is (N, M)).
      image : 2D array, dtype=float, optional
          Image containing an initial reconstruction estimate. Shape of this
          array should be ``(radon_image.shape[0], radon_image.shape[0])``. The
          default is an array of zeros.
      projection_shifts : 1D array, dtype=float
          Shift the projections contained in ``radon_image`` (the sinogram) by
          this many pixels before reconstructing the image. The i'th value
          defines the shift of the i'th column of ``radon_image``.
      clip : length-2 sequence of floats
          Force all values in the reconstructed tomogram to lie in the range
          ``[clip[0], clip[1]]``
      relaxation : float
          Relaxation parameter for the update step. A higher value can
          improve the convergence rate, but one runs the risk of instabilities.
          Values close to or higher than 1 are not recommended.
  
      Returns
      -------
      reconstructed : ndarray
          Reconstructed image. The rotation axis will be located in the pixel
          with indices
          ``(reconstructed.shape[0] // 2, reconstructed.shape[1] // 2)``.
  
      Notes
      -----
      Algebraic Reconstruction Techniques are based on formulating the tomography
      reconstruction problem as a set of linear equations. Along each ray,
      the projected value is the sum of all the values of the cross section along
      the ray. A typical feature of SART (and a few other variants of algebraic
      techniques) is that it samples the cross section at equidistant points
      along the ray, using linear interpolation between the pixel values of the
      cross section. The resulting set of linear equations are then solved using
      a slightly modified Kaczmarz method.
  
      When using SART, a single iteration is usually sufficient to obtain a good
      reconstruction. Further iterations will tend to enhance high-frequency
      information, but will also often increase the noise.
  
      References
      ----------
      .. [1] AC Kak, M Slaney, "Principles of Computerized Tomographic
             Imaging", IEEE Press 1988.
      .. [2] AH Andersen, AC Kak, "Simultaneous algebraic reconstruction
             technique (SART): a superior implementation of the ART algorithm",
             Ultrasonic Imaging 6 pp 81--94 (1984)
      .. [3] S Kaczmarz, "Angenäherte auflösung von systemen linearer
             gleichungen", Bulletin International de l’Academie Polonaise des
             Sciences et des Lettres 35 pp 355--357 (1937)
      .. [4] Kohler, T. "A projection access scheme for iterative
             reconstruction based on the golden section." Nuclear Science
             Symposium Conference Record, 2004 IEEE. Vol. 6. IEEE, 2004.
      .. [5] Kaczmarz' method, Wikipedia,
             http://en.wikipedia.org/wiki/Kaczmarz_method
      """
      if radon_image.ndim != 2:
          raise ValueError('radon_image must be two dimensional')
      reconstructed_shape = (radon_image.shape[0], radon_image.shape[0])
      if theta is None:
          theta = np.linspace(0, 180, radon_image.shape[1], endpoint=False)
      elif theta.shape != (radon_image.shape[1],):
          raise ValueError('Shape of theta (%s) does not match the '
                           'number of projections (%d)'
                           % (projection_shifts.shape, radon_image.shape[1]))
      if image is None:
          image = np.zeros(reconstructed_shape, dtype=np.float)
      elif image.shape != reconstructed_shape:
          raise ValueError('Shape of image (%s) does not match first dimension '
                           'of radon_image (%s)'
                           % (image.shape, reconstructed_shape))
      if projection_shifts is None:
          projection_shifts = np.zeros((radon_image.shape[1],), dtype=np.float)
      elif projection_shifts.shape != (radon_image.shape[1],):
          raise ValueError('Shape of projection_shifts (%s) does not match the '
                           'number of projections (%d)'
                           % (projection_shifts.shape, radon_image.shape[1]))
      if not clip is None:
          if len(clip) != 2:
              raise ValueError('clip must be a length-2 sequence')
          clip = (float(clip[0]), float(clip[1]))
      relaxation = float(relaxation)
  
      for angle_index in order_angles_golden_ratio(theta):
          image_update = sart_projection_update(image, theta[angle_index],
                                                radon_image[:, angle_index],
                                                projection_shifts[angle_index])
          image += relaxation * image_update
          if not clip is None:
              image = np.clip(image, clip[0], clip[1])
      return image
   from numpy.testing import assert_array_equal, assert_almost_equal, \
      assert_array_almost_equal, assert_raises, assert_equal
  import numpy as np
  import math
  
  from skimage.measure._regionprops import (regionprops, PROPS, perimeter,
                                            _parse_docs)
  
  
  SAMPLE = np.array(
      [[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 1, 1, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
       [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
       [1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, 1, 0],
       [0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, 1, 1],
       [0, 1, 1, 0, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1]]
  )
  INTENSITY_SAMPLE = SAMPLE.copy()
  INTENSITY_SAMPLE[1, 9:11] = 2
  
  SAMPLE_3D = np.zeros((6, 6, 6), dtype=np.uint8)
  SAMPLE_3D[1:3, 1:3, 1:3] = 1
  SAMPLE_3D[3, 2, 2] = 1 
  INTENSITY_SAMPLE_3D = SAMPLE_3D.copy()
  
  def test_all_props():
      region = regionprops(SAMPLE, INTENSITY_SAMPLE)[0]
      for prop in PROPS:
          assert_almost_equal(region[prop], getattr(region, PROPS[prop]))
  
  
  def test_all_props_3d():
      region = regionprops(SAMPLE_3D, INTENSITY_SAMPLE_3D)[0]
      for prop in PROPS:
          try:
              assert_almost_equal(region[prop], getattr(region, PROPS[prop]))
          except NotImplementedError:
              pass
  
  def test_dtype():
      regionprops(np.zeros((10, 10), dtype=np.int))
      regionprops(np.zeros((10, 10), dtype=np.uint))
      assert_raises((TypeError), regionprops,
                    np.zeros((10, 10), dtype=np.float))
      assert_raises((TypeError), regionprops,
                    np.zeros((10, 10), dtype=np.double))
  
  
  def test_ndim():
      regionprops(np.zeros((10, 10), dtype=np.int))
      regionprops(np.zeros((10, 10, 1), dtype=np.int))
      regionprops(np.zeros((10, 10, 1, 1), dtype=np.int))
      regionprops(np.zeros((10, 10, 10), dtype=np.int))
      assert_raises(TypeError, regionprops, np.zeros((10, 10, 10, 2), dtype=np.int))
  
  
  def test_area():
      area = regionprops(SAMPLE)[0].area
      assert area == np.sum(SAMPLE)
      area = regionprops(SAMPLE_3D)[0].area
      assert area == np.sum(SAMPLE_3D)
  
  
  def test_bbox():
      bbox = regionprops(SAMPLE)[0].bbox
      assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1]))
  
      SAMPLE_mod = SAMPLE.copy()
      SAMPLE_mod[:, -1] = 0
      bbox = regionprops(SAMPLE_mod)[0].bbox
      assert_array_almost_equal(bbox, (0, 0, SAMPLE.shape[0], SAMPLE.shape[1]-1))
  
      bbox = regionprops(SAMPLE_3D)[0].bbox
      assert_array_almost_equal(bbox, (1, 1, 1, 4, 3, 3))
  
  
  def test_moments_central():
      mu = regionprops(SAMPLE)[0].moments_central
      # determined with OpenCV
      assert_almost_equal(mu[0,2], 436.00000000000045)
      # different from OpenCV results, bug in OpenCV
      assert_almost_equal(mu[0,3], -737.333333333333)
      assert_almost_equal(mu[1,1], -87.33333333333303)
      assert_almost_equal(mu[1,2], -127.5555555555593)
      assert_almost_equal(mu[2,0], 1259.7777777777774)
      assert_almost_equal(mu[2,1], 2000.296296296291)
      assert_almost_equal(mu[3,0], -760.0246913580195)
  
  
  def test_centroid():
      centroid = regionprops(SAMPLE)[0].centroid
      # determined with MATLAB
      assert_array_almost_equal(centroid, (5.66666666666666, 9.444444444444444))
  
  
  def test_convex_area():
      area = regionprops(SAMPLE)[0].convex_area
      # determined with MATLAB
      assert area == 124
  
  
  def test_convex_image():
      img = regionprops(SAMPLE)[0].convex_image
      # determined with MATLAB
      ref = np.array(
          [[0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0],
           [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0],
           [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
           [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0],
           [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0],
           [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
           [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0],
           [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]]
      )
      assert_array_equal(img, ref)
  
  
  def test_coordinates():
      sample = np.zeros((10, 10), dtype=np.int8)
      coords = np.array([[3, 2], [3, 3], [3, 4]])
      sample[coords[:, 0], coords[:, 1]] = 1
      prop_coords = regionprops(sample)[0].coords
      assert_array_equal(prop_coords, coords)
  
      sample = np.zeros((6, 6, 6), dtype=np.int8)
      coords = np.array([[1, 1, 1], [1, 2, 1], [1, 3, 1]])
      sample[coords[:, 0], coords[:, 1], coords[:, 2]] = 1
      prop_coords = regionprops(sample)[0].coords
      assert_array_equal(prop_coords, coords)
  
  def test_eccentricity():
      eps = regionprops(SAMPLE)[0].eccentricity
      assert_almost_equal(eps, 0.814629313427)
  
      img = np.zeros((5, 5), dtype=np.int)
      img[2, 2] = 1
      eps = regionprops(img)[0].eccentricity
      assert_almost_equal(eps, 0)
  
  
  def test_equiv_diameter():
      diameter = regionprops(SAMPLE)[0].equivalent_diameter
      # determined with MATLAB
      assert_almost_equal(diameter, 9.57461472963)
  
  
  def test_euler_number():
      en = regionprops(SAMPLE)[0].euler_number
      assert en == 1
  
      SAMPLE_mod = SAMPLE.copy()
      SAMPLE_mod[7, -3] = 0
      en = regionprops(SAMPLE_mod)[0].euler_number
      assert en == 0
  
  
  def test_extent():
      extent = regionprops(SAMPLE)[0].extent
      assert_almost_equal(extent, 0.4)
  
  
  def test_moments_hu():
      hu = regionprops(SAMPLE)[0].moments_hu
      ref = np.array([
           3.27117627e-01,
           2.63869194e-02,
           2.35390060e-02,
           1.23151193e-03,
           1.38882330e-06,
          -2.72586158e-05,
           6.48350653e-06
      ])
      # bug in OpenCV caused in Central Moments calculation?
      assert_array_almost_equal(hu, ref)
  
  
  def test_image():
      img = regionprops(SAMPLE)[0].image
      assert_array_equal(img, SAMPLE)
  
      img = regionprops(SAMPLE_3D)[0].image
      assert_array_equal(img, SAMPLE_3D[1:4, 1:3, 1:3])
  
  
  def test_label():
      label = regionprops(SAMPLE)[0].label
      assert_array_equal(label, 1)
  
      label = regionprops(SAMPLE_3D)[0].label
      assert_array_equal(label, 1)
  
  
  def test_filled_area():
      area = regionprops(SAMPLE)[0].filled_area
      assert area == np.sum(SAMPLE)
  
      SAMPLE_mod = SAMPLE.copy()
      SAMPLE_mod[7, -3] = 0
      area = regionprops(SAMPLE_mod)[0].filled_area
      assert area == np.sum(SAMPLE)
  
  
  def test_filled_image():
      img = regionprops(SAMPLE)[0].filled_image
      assert_array_equal(img, SAMPLE)
  
  
  def test_major_axis_length():
      length = regionprops(SAMPLE)[0].major_axis_length
      # MATLAB has different interpretation of ellipse than found in literature,
      # here implemented as found in literature
      assert_almost_equal(length, 16.7924234999)
  
  
  def test_max_intensity():
      intensity = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE
                              )[0].max_intensity
      assert_almost_equal(intensity, 2)
  
  
  def test_mean_intensity():
      intensity = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE
                              )[0].mean_intensity
      assert_almost_equal(intensity, 1.02777777777777)
  
  
  def test_min_intensity():
      intensity = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE
                              )[0].min_intensity
      assert_almost_equal(intensity, 1)
  
  
  def test_minor_axis_length():
      length = regionprops(SAMPLE)[0].minor_axis_length
      # MATLAB has different interpretation of ellipse than found in literature,
      # here implemented as found in literature
      assert_almost_equal(length, 9.739302807263)
  
  
  def test_moments():
      m = regionprops(SAMPLE)[0].moments
      # determined with OpenCV
      assert_almost_equal(m[0,0], 72.0)
      assert_almost_equal(m[0,1], 408.0)
      assert_almost_equal(m[0,2], 2748.0)
      assert_almost_equal(m[0,3], 19776.0)
      assert_almost_equal(m[1,0], 680.0)
      assert_almost_equal(m[1,1], 3766.0)
      assert_almost_equal(m[1,2], 24836.0)
      assert_almost_equal(m[2,0], 7682.0)
      assert_almost_equal(m[2,1], 43882.0)
      assert_almost_equal(m[3,0], 95588.0)
  
  
  def test_moments_normalized():
      nu = regionprops(SAMPLE)[0].moments_normalized
      # determined with OpenCV
      assert_almost_equal(nu[0,2], 0.08410493827160502)
      assert_almost_equal(nu[1,1], -0.016846707818929982)
      assert_almost_equal(nu[1,2], -0.002899800614433943)
      assert_almost_equal(nu[2,0], 0.24301268861454037)
      assert_almost_equal(nu[2,1], 0.045473992910668816)
      assert_almost_equal(nu[3,0], -0.017278118992041805)
  
  
  def test_orientation():
      orientation = regionprops(SAMPLE)[0].orientation
      # determined with MATLAB
      assert_almost_equal(orientation, 0.10446844651921)
      # test correct quadrant determination
      orientation2 = regionprops(SAMPLE.T)[0].orientation
      assert_almost_equal(orientation2, math.pi / 2 - orientation)
      # test diagonal regions
      diag = np.eye(10, dtype=int)
      orientation_diag = regionprops(diag)[0].orientation
      assert_almost_equal(orientation_diag, -math.pi / 4)
      orientation_diag = regionprops(np.flipud(diag))[0].orientation
      assert_almost_equal(orientation_diag, math.pi / 4)
      orientation_diag = regionprops(np.fliplr(diag))[0].orientation
      assert_almost_equal(orientation_diag, math.pi / 4)
      orientation_diag = regionprops(np.fliplr(np.flipud(diag)))[0].orientation
      assert_almost_equal(orientation_diag, -math.pi / 4)
  
  
  def test_perimeter():
      per = regionprops(SAMPLE)[0].perimeter
      assert_almost_equal(per, 55.2487373415)
  
      per = perimeter(SAMPLE.astype('double'), neighbourhood=8)
      assert_almost_equal(per, 46.8284271247)
  
  
  def test_solidity():
      solidity = regionprops(SAMPLE)[0].solidity
      # determined with MATLAB
      assert_almost_equal(solidity, 0.580645161290323)
  
  
  def test_weighted_moments_central():
      wmu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE
                        )[0].weighted_moments_central
      ref = np.array(
          [[  7.4000000000e+01, -2.1316282073e-13,  4.7837837838e+02,
              -7.5943608473e+02],
           [  3.7303493627e-14, -8.7837837838e+01, -1.4801314828e+02,
              -1.2714707125e+03],
           [  1.2602837838e+03,  2.1571526662e+03,  6.6989799420e+03,
               1.5304076361e+04],
           [ -7.6561796932e+02, -4.2385971907e+03, -9.9501164076e+03,
              -3.3156729271e+04]]
      )
      np.set_printoptions(precision=10)
      assert_array_almost_equal(wmu, ref)
  
  
  def test_weighted_centroid():
      centroid = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE
                             )[0].weighted_centroid
      assert_array_almost_equal(centroid, (5.540540540540, 9.445945945945))
  
  
  def test_weighted_moments_hu():
      whu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE
                        )[0].weighted_moments_hu
      ref = np.array([
          3.1750587329e-01,
          2.1417517159e-02,
          2.3609322038e-02,
          1.2565683360e-03,
          8.3014209421e-07,
          -3.5073773473e-05,
          6.7936409056e-06
      ])
      assert_array_almost_equal(whu, ref)
  
  
  def test_weighted_moments():
      wm = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE
                       )[0].weighted_moments
      ref = np.array(
          [[  7.4000000000e+01, 4.1000000000e+02, 2.7500000000e+03,
              1.9778000000e+04],
           [  6.9900000000e+02, 3.7850000000e+03, 2.4855000000e+04,
              1.7500100000e+05],
           [  7.8630000000e+03, 4.4063000000e+04, 2.9347700000e+05,
              2.0810510000e+06],
           [  9.7317000000e+04, 5.7256700000e+05, 3.9007170000e+06,
              2.8078871000e+07]]
      )
      assert_array_almost_equal(wm, ref)
  
  
  def test_weighted_moments_normalized():
      wnu = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE
                        )[0].weighted_moments_normalized
      ref = np.array(
          [[       np.nan,        np.nan,  0.0873590903, -0.0161217406],
           [       np.nan, -0.0160405109, -0.0031421072, -0.0031376984],
           [  0.230146783,  0.0457932622,  0.0165315478,  0.0043903193],
           [-0.0162529732, -0.0104598869, -0.0028544152, -0.0011057191]]
      )
      assert_array_almost_equal(wnu, ref)
  
  
  def test_label_sequence():
      a = np.empty((2, 2), dtype=np.int)
      a[:, :] = 2
      ps = regionprops(a)
      assert len(ps) == 1
      assert ps[0].label == 2
  
  
  def test_pure_background():
      a = np.zeros((2, 2), dtype=np.int)
      ps = regionprops(a)
      assert len(ps) == 0
  
  
  def test_invalid():
      ps = regionprops(SAMPLE)
  
      def get_intensity_image():
          ps[0].intensity_image
  
      assert_raises(AttributeError, get_intensity_image)
  
  
  def test_invalid_size():
      wrong_intensity_sample = np.array([[1], [1]])
      assert_raises(ValueError, regionprops, SAMPLE, wrong_intensity_sample)
  
  
  def test_equals():
      arr = np.zeros((100, 100), dtype=np.int)
      arr[0:25, 0:25] = 1
      arr[50:99, 50:99] = 2
  
      regions = regionprops(arr)
      r1 = regions[0]
  
      regions = regionprops(arr)
      r2 = regions[0]
      r3 = regions[1]
  
      assert_equal(r1 == r2, True, "Same regionprops are not equal")
      assert_equal(r1 != r3, True, "Different regionprops are equal")
  
  
  def test_iterate_all_props():
      region = regionprops(SAMPLE)[0]
      p0 = dict((p, region[p]) for p in region)
  
      region = regionprops(SAMPLE, intensity_image=INTENSITY_SAMPLE)[0]
      p1 = dict((p, region[p]) for p in region)
  
      assert len(p0) < len(p1)
  
  
  def test_cache():
      region = regionprops(SAMPLE)[0]
      f0 = region.filled_image
      region._label_image[:10] = 1
      f1 = region.filled_image
  
      # Changed underlying image, but cache keeps result the same
      assert_array_equal(f0, f1)
  
      # Now invalidate cache
      region._cache_active = False
      f1 = region.filled_image
  
      assert np.any(f0 != f1)
  
  
  def test_docstrings_and_props():
      region = regionprops(SAMPLE)[0]
  
      docs = _parse_docs()
      props = [m for m in dir(region) if not m.startswith('_')]
  
      nr_docs_parsed = len(docs)
      nr_props = len(props)
      assert_equal(nr_docs_parsed, nr_props)
  
      ds = docs['weighted_moments_normalized']
      assert 'iteration' not in ds
      assert len(ds.split('\n')) > 3
  
  
  if __name__ == "__main__":
      from numpy.testing import run_module_suite
      run_module_suite()
  
apis: scipy.interpolate.interp1d numpy.argsort numpy.asarray numpy.abs numpy.cos numpy.linspace numpy.deg2rad scipy.fftpack.ifft numpy.ceil numpy.zeros numpy.searchsorted numpy.arange numpy.all numpy.array numpy.log2 numpy.interp scipy.fftpack.fft scipy.fftpack.fftfreq numpy.clip numpy.sqrt numpy.sin
apis: numpy.testing.assert_almost_equal numpy.eye numpy.testing.assert_raises numpy.sum numpy.empty numpy.zeros numpy.flipud numpy.testing.assert_equal numpy.testing.run_module_suite numpy.any numpy.set_printoptions numpy.testing.assert_array_equal numpy.fliplr numpy.testing.assert_array_almost_equal numpy.array