| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365 |
- import numpy as np
- from scipy.ndimage import distance_transform_edt as distance
- from .._shared.utils import _supported_float_type
- def _cv_calculate_variation(image, phi, mu, lambda1, lambda2, dt):
- """Returns the variation of level set 'phi' based on algorithm parameters.
- This corresponds to equation (22) of the paper by Pascal Getreuer,
- which computes the next iteration of the level set based on a current
- level set.
- A full explanation regarding all the terms is beyond the scope of the
- present description, but there is one difference of particular import.
- In the original algorithm, convergence is accelerated, and required
- memory is reduced, by using a single array. This array, therefore, is a
- combination of non-updated and updated values. If this were to be
- implemented in python, this would require a double loop, where the
- benefits of having fewer iterations would be outweided by massively
- increasing the time required to perform each individual iteration. A
- similar approach is used by Rami Cohen, and it is from there that the
- C1-4 notation is taken.
- """
- eta = 1e-16
- P = np.pad(phi, 1, mode='edge')
- phixp = P[1:-1, 2:] - P[1:-1, 1:-1]
- phixn = P[1:-1, 1:-1] - P[1:-1, :-2]
- phix0 = (P[1:-1, 2:] - P[1:-1, :-2]) / 2.0
- phiyp = P[2:, 1:-1] - P[1:-1, 1:-1]
- phiyn = P[1:-1, 1:-1] - P[:-2, 1:-1]
- phiy0 = (P[2:, 1:-1] - P[:-2, 1:-1]) / 2.0
- C1 = 1.0 / np.sqrt(eta + phixp**2 + phiy0**2)
- C2 = 1.0 / np.sqrt(eta + phixn**2 + phiy0**2)
- C3 = 1.0 / np.sqrt(eta + phix0**2 + phiyp**2)
- C4 = 1.0 / np.sqrt(eta + phix0**2 + phiyn**2)
- K = P[1:-1, 2:] * C1 + P[1:-1, :-2] * C2 + P[2:, 1:-1] * C3 + P[:-2, 1:-1] * C4
- Hphi = (phi > 0).astype(image.dtype)
- (c1, c2) = _cv_calculate_averages(image, Hphi)
- difference_from_average_term = (
- -lambda1 * (image - c1) ** 2 + lambda2 * (image - c2) ** 2
- )
- new_phi = phi + (dt * _cv_delta(phi)) * (mu * K + difference_from_average_term)
- return new_phi / (1 + mu * dt * _cv_delta(phi) * (C1 + C2 + C3 + C4))
- def _cv_heavyside(x, eps=1.0):
- """Returns the result of a regularised heavyside function of the
- input value(s).
- """
- return 0.5 * (1.0 + (2.0 / np.pi) * np.arctan(x / eps))
- def _cv_delta(x, eps=1.0):
- """Returns the result of a regularised dirac function of the
- input value(s).
- """
- return eps / (eps**2 + x**2)
- def _cv_calculate_averages(image, Hphi):
- """Returns the average values 'inside' and 'outside'."""
- H = Hphi
- Hinv = 1.0 - H
- Hsum = np.sum(H)
- Hinvsum = np.sum(Hinv)
- avg_inside = np.sum(image * H)
- avg_oustide = np.sum(image * Hinv)
- if Hsum != 0:
- avg_inside /= Hsum
- if Hinvsum != 0:
- avg_oustide /= Hinvsum
- return (avg_inside, avg_oustide)
- def _cv_difference_from_average_term(image, Hphi, lambda_pos, lambda_neg):
- """Returns the 'energy' contribution due to the difference from
- the average value within a region at each point.
- """
- (c1, c2) = _cv_calculate_averages(image, Hphi)
- Hinv = 1.0 - Hphi
- return lambda_pos * (image - c1) ** 2 * Hphi + lambda_neg * (image - c2) ** 2 * Hinv
- def _cv_edge_length_term(phi, mu):
- """Returns the 'energy' contribution due to the length of the
- edge between regions at each point, multiplied by a factor 'mu'.
- """
- P = np.pad(phi, 1, mode='edge')
- fy = (P[2:, 1:-1] - P[:-2, 1:-1]) / 2.0
- fx = (P[1:-1, 2:] - P[1:-1, :-2]) / 2.0
- return mu * _cv_delta(phi) * np.sqrt(fx**2 + fy**2)
- def _cv_energy(image, phi, mu, lambda1, lambda2):
- """Returns the total 'energy' of the current level set function.
- This corresponds to equation (7) of the paper by Pascal Getreuer,
- which is the weighted sum of the following:
- (A) the length of the contour produced by the zero values of the
- level set,
- (B) the area of the "foreground" (area of the image where the
- level set is positive),
- (C) the variance of the image inside the foreground,
- (D) the variance of the image outside of the foreground
- Each value is computed for each pixel, and then summed. The weight
- of (B) is set to 0 in this implementation.
- """
- H = _cv_heavyside(phi)
- avgenergy = _cv_difference_from_average_term(image, H, lambda1, lambda2)
- lenenergy = _cv_edge_length_term(phi, mu)
- return np.sum(avgenergy) + np.sum(lenenergy)
- def _cv_reset_level_set(phi):
- """This is a placeholder function as resetting the level set is not
- strictly necessary, and has not been done for this implementation.
- """
- return phi
- def _cv_checkerboard(image_size, square_size, dtype=np.float64):
- """Generates a checkerboard level set function.
- According to Pascal Getreuer, such a level set function has fast
- convergence.
- """
- yv = np.arange(image_size[0], dtype=dtype).reshape(image_size[0], 1)
- xv = np.arange(image_size[1], dtype=dtype)
- sf = np.pi / square_size
- xv *= sf
- yv *= sf
- return np.sin(yv) * np.sin(xv)
- def _cv_large_disk(image_size):
- """Generates a disk level set function.
- The disk covers the whole image along its smallest dimension.
- """
- res = np.ones(image_size)
- centerY = int((image_size[0] - 1) / 2)
- centerX = int((image_size[1] - 1) / 2)
- res[centerY, centerX] = 0.0
- radius = float(min(centerX, centerY))
- return (radius - distance(res)) / radius
- def _cv_small_disk(image_size):
- """Generates a disk level set function.
- The disk covers half of the image along its smallest dimension.
- """
- res = np.ones(image_size)
- centerY = int((image_size[0] - 1) / 2)
- centerX = int((image_size[1] - 1) / 2)
- res[centerY, centerX] = 0.0
- radius = float(min(centerX, centerY)) / 2.0
- return (radius - distance(res)) / (radius * 3)
- def _cv_init_level_set(init_level_set, image_shape, dtype=np.float64):
- """Generates an initial level set function conditional on input arguments."""
- if isinstance(init_level_set, str):
- if init_level_set == 'checkerboard':
- res = _cv_checkerboard(image_shape, 5, dtype)
- elif init_level_set == 'disk':
- res = _cv_large_disk(image_shape)
- elif init_level_set == 'small disk':
- res = _cv_small_disk(image_shape)
- else:
- raise ValueError("Incorrect name for starting level set preset.")
- else:
- res = init_level_set
- return res.astype(dtype, copy=False)
- def chan_vese(
- image,
- mu=0.25,
- lambda1=1.0,
- lambda2=1.0,
- tol=1e-3,
- max_num_iter=500,
- dt=0.5,
- init_level_set='checkerboard',
- extended_output=False,
- ):
- """Chan-Vese segmentation algorithm.
- Active contour model by evolving a level set. Can be used to
- segment objects without clearly defined boundaries.
- Parameters
- ----------
- image : (M, N) ndarray
- Grayscale image to be segmented.
- mu : float, optional
- 'edge length' weight parameter. Higher `mu` values will
- produce a 'round' edge, while values closer to zero will
- detect smaller objects.
- lambda1 : float, optional
- 'difference from average' weight parameter for the output
- region with value 'True'. If it is lower than `lambda2`, this
- region will have a larger range of values than the other.
- lambda2 : float, optional
- 'difference from average' weight parameter for the output
- region with value 'False'. If it is lower than `lambda1`, this
- region will have a larger range of values than the other.
- tol : float, positive, optional
- Level set variation tolerance between iterations. If the
- L2 norm difference between the level sets of successive
- iterations normalized by the area of the image is below this
- value, the algorithm will assume that the solution was
- reached.
- max_num_iter : uint, optional
- Maximum number of iterations allowed before the algorithm
- interrupts itself.
- dt : float, optional
- A multiplication factor applied at calculations for each step,
- serves to accelerate the algorithm. While higher values may
- speed up the algorithm, they may also lead to convergence
- problems.
- init_level_set : str or (M, N) ndarray, optional
- Defines the starting level set used by the algorithm.
- If a string is inputted, a level set that matches the image
- size will automatically be generated. Alternatively, it is
- possible to define a custom level set, which should be an
- array of float values, with the same shape as 'image'.
- Accepted string values are as follows.
- 'checkerboard'
- the starting level set is defined as
- sin(x/5*pi)*sin(y/5*pi), where x and y are pixel
- coordinates. This level set has fast convergence, but may
- fail to detect implicit edges.
- 'disk'
- the starting level set is defined as the opposite
- of the distance from the center of the image minus half of
- the minimum value between image width and image height.
- This is somewhat slower, but is more likely to properly
- detect implicit edges.
- 'small disk'
- the starting level set is defined as the
- opposite of the distance from the center of the image
- minus a quarter of the minimum value between image width
- and image height.
- extended_output : bool, optional
- If set to True, the return value will be a tuple containing
- the three return values (see below). If set to False which
- is the default value, only the 'segmentation' array will be
- returned.
- Returns
- -------
- segmentation : (M, N) ndarray, bool
- Segmentation produced by the algorithm.
- phi : (M, N) ndarray of floats
- Final level set computed by the algorithm.
- energies : list of floats
- Shows the evolution of the 'energy' for each step of the
- algorithm. This should allow to check whether the algorithm
- converged.
- Notes
- -----
- The Chan-Vese Algorithm is designed to segment objects without
- clearly defined boundaries. This algorithm is based on level sets
- that are evolved iteratively to minimize an energy, which is
- defined by weighted values corresponding to the sum of differences
- intensity from the average value outside the segmented region, the
- sum of differences from the average value inside the segmented
- region, and a term which is dependent on the length of the
- boundary of the segmented region.
- This algorithm was first proposed by Tony Chan and Luminita Vese,
- in a publication entitled "An Active Contour Model Without Edges"
- [1]_.
- This implementation of the algorithm is somewhat simplified in the
- sense that the area factor 'nu' described in the original paper is
- not implemented, and is only suitable for grayscale images.
- Typical values for `lambda1` and `lambda2` are 1. If the
- 'background' is very different from the segmented object in terms
- of distribution (for example, a uniform black image with figures
- of varying intensity), then these values should be different from
- each other.
- Typical values for mu are between 0 and 1, though higher values
- can be used when dealing with shapes with very ill-defined
- contours.
- The 'energy' which this algorithm tries to minimize is defined
- as the sum of the differences from the average within the region
- squared and weighed by the 'lambda' factors to which is added the
- length of the contour multiplied by the 'mu' factor.
- Supports 2D grayscale images only, and does not implement the area
- term described in the original article.
- References
- ----------
- .. [1] An Active Contour Model without Edges, Tony Chan and
- Luminita Vese, Scale-Space Theories in Computer Vision,
- 1999, :DOI:`10.1007/3-540-48236-9_13`
- .. [2] Chan-Vese Segmentation, Pascal Getreuer Image Processing On
- Line, 2 (2012), pp. 214-224,
- :DOI:`10.5201/ipol.2012.g-cv`
- .. [3] The Chan-Vese Algorithm - Project Report, Rami Cohen, 2011
- :arXiv:`1107.2782`
- """
- if len(image.shape) != 2:
- raise ValueError("Input image should be a 2D array.")
- float_dtype = _supported_float_type(image.dtype)
- phi = _cv_init_level_set(init_level_set, image.shape, dtype=float_dtype)
- if type(phi) != np.ndarray or phi.shape != image.shape:
- raise ValueError(
- "The dimensions of initial level set do not "
- "match the dimensions of image."
- )
- image = image.astype(float_dtype, copy=False)
- image = image - np.min(image)
- if np.max(image) != 0:
- image = image / np.max(image)
- i = 0
- old_energy = _cv_energy(image, phi, mu, lambda1, lambda2)
- energies = []
- phivar = tol + 1
- segmentation = phi > 0
- while phivar > tol and i < max_num_iter:
- # Save old level set values
- oldphi = phi
- # Calculate new level set
- phi = _cv_calculate_variation(image, phi, mu, lambda1, lambda2, dt)
- phi = _cv_reset_level_set(phi)
- phivar = np.sqrt(((phi - oldphi) ** 2).mean())
- # Extract energy and compare to previous level set and
- # segmentation to see if continuing is necessary
- segmentation = phi > 0
- new_energy = _cv_energy(image, phi, mu, lambda1, lambda2)
- # Save old energy values
- energies.append(old_energy)
- old_energy = new_energy
- i += 1
- if extended_output:
- return (segmentation, phi, energies)
- else:
- return segmentation
|