_chan_vese.py 13 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365
  1. import numpy as np
  2. from scipy.ndimage import distance_transform_edt as distance
  3. from .._shared.utils import _supported_float_type
  4. def _cv_calculate_variation(image, phi, mu, lambda1, lambda2, dt):
  5. """Returns the variation of level set 'phi' based on algorithm parameters.
  6. This corresponds to equation (22) of the paper by Pascal Getreuer,
  7. which computes the next iteration of the level set based on a current
  8. level set.
  9. A full explanation regarding all the terms is beyond the scope of the
  10. present description, but there is one difference of particular import.
  11. In the original algorithm, convergence is accelerated, and required
  12. memory is reduced, by using a single array. This array, therefore, is a
  13. combination of non-updated and updated values. If this were to be
  14. implemented in python, this would require a double loop, where the
  15. benefits of having fewer iterations would be outweided by massively
  16. increasing the time required to perform each individual iteration. A
  17. similar approach is used by Rami Cohen, and it is from there that the
  18. C1-4 notation is taken.
  19. """
  20. eta = 1e-16
  21. P = np.pad(phi, 1, mode='edge')
  22. phixp = P[1:-1, 2:] - P[1:-1, 1:-1]
  23. phixn = P[1:-1, 1:-1] - P[1:-1, :-2]
  24. phix0 = (P[1:-1, 2:] - P[1:-1, :-2]) / 2.0
  25. phiyp = P[2:, 1:-1] - P[1:-1, 1:-1]
  26. phiyn = P[1:-1, 1:-1] - P[:-2, 1:-1]
  27. phiy0 = (P[2:, 1:-1] - P[:-2, 1:-1]) / 2.0
  28. C1 = 1.0 / np.sqrt(eta + phixp**2 + phiy0**2)
  29. C2 = 1.0 / np.sqrt(eta + phixn**2 + phiy0**2)
  30. C3 = 1.0 / np.sqrt(eta + phix0**2 + phiyp**2)
  31. C4 = 1.0 / np.sqrt(eta + phix0**2 + phiyn**2)
  32. K = P[1:-1, 2:] * C1 + P[1:-1, :-2] * C2 + P[2:, 1:-1] * C3 + P[:-2, 1:-1] * C4
  33. Hphi = (phi > 0).astype(image.dtype)
  34. (c1, c2) = _cv_calculate_averages(image, Hphi)
  35. difference_from_average_term = (
  36. -lambda1 * (image - c1) ** 2 + lambda2 * (image - c2) ** 2
  37. )
  38. new_phi = phi + (dt * _cv_delta(phi)) * (mu * K + difference_from_average_term)
  39. return new_phi / (1 + mu * dt * _cv_delta(phi) * (C1 + C2 + C3 + C4))
  40. def _cv_heavyside(x, eps=1.0):
  41. """Returns the result of a regularised heavyside function of the
  42. input value(s).
  43. """
  44. return 0.5 * (1.0 + (2.0 / np.pi) * np.arctan(x / eps))
  45. def _cv_delta(x, eps=1.0):
  46. """Returns the result of a regularised dirac function of the
  47. input value(s).
  48. """
  49. return eps / (eps**2 + x**2)
  50. def _cv_calculate_averages(image, Hphi):
  51. """Returns the average values 'inside' and 'outside'."""
  52. H = Hphi
  53. Hinv = 1.0 - H
  54. Hsum = np.sum(H)
  55. Hinvsum = np.sum(Hinv)
  56. avg_inside = np.sum(image * H)
  57. avg_oustide = np.sum(image * Hinv)
  58. if Hsum != 0:
  59. avg_inside /= Hsum
  60. if Hinvsum != 0:
  61. avg_oustide /= Hinvsum
  62. return (avg_inside, avg_oustide)
  63. def _cv_difference_from_average_term(image, Hphi, lambda_pos, lambda_neg):
  64. """Returns the 'energy' contribution due to the difference from
  65. the average value within a region at each point.
  66. """
  67. (c1, c2) = _cv_calculate_averages(image, Hphi)
  68. Hinv = 1.0 - Hphi
  69. return lambda_pos * (image - c1) ** 2 * Hphi + lambda_neg * (image - c2) ** 2 * Hinv
  70. def _cv_edge_length_term(phi, mu):
  71. """Returns the 'energy' contribution due to the length of the
  72. edge between regions at each point, multiplied by a factor 'mu'.
  73. """
  74. P = np.pad(phi, 1, mode='edge')
  75. fy = (P[2:, 1:-1] - P[:-2, 1:-1]) / 2.0
  76. fx = (P[1:-1, 2:] - P[1:-1, :-2]) / 2.0
  77. return mu * _cv_delta(phi) * np.sqrt(fx**2 + fy**2)
  78. def _cv_energy(image, phi, mu, lambda1, lambda2):
  79. """Returns the total 'energy' of the current level set function.
  80. This corresponds to equation (7) of the paper by Pascal Getreuer,
  81. which is the weighted sum of the following:
  82. (A) the length of the contour produced by the zero values of the
  83. level set,
  84. (B) the area of the "foreground" (area of the image where the
  85. level set is positive),
  86. (C) the variance of the image inside the foreground,
  87. (D) the variance of the image outside of the foreground
  88. Each value is computed for each pixel, and then summed. The weight
  89. of (B) is set to 0 in this implementation.
  90. """
  91. H = _cv_heavyside(phi)
  92. avgenergy = _cv_difference_from_average_term(image, H, lambda1, lambda2)
  93. lenenergy = _cv_edge_length_term(phi, mu)
  94. return np.sum(avgenergy) + np.sum(lenenergy)
  95. def _cv_reset_level_set(phi):
  96. """This is a placeholder function as resetting the level set is not
  97. strictly necessary, and has not been done for this implementation.
  98. """
  99. return phi
  100. def _cv_checkerboard(image_size, square_size, dtype=np.float64):
  101. """Generates a checkerboard level set function.
  102. According to Pascal Getreuer, such a level set function has fast
  103. convergence.
  104. """
  105. yv = np.arange(image_size[0], dtype=dtype).reshape(image_size[0], 1)
  106. xv = np.arange(image_size[1], dtype=dtype)
  107. sf = np.pi / square_size
  108. xv *= sf
  109. yv *= sf
  110. return np.sin(yv) * np.sin(xv)
  111. def _cv_large_disk(image_size):
  112. """Generates a disk level set function.
  113. The disk covers the whole image along its smallest dimension.
  114. """
  115. res = np.ones(image_size)
  116. centerY = int((image_size[0] - 1) / 2)
  117. centerX = int((image_size[1] - 1) / 2)
  118. res[centerY, centerX] = 0.0
  119. radius = float(min(centerX, centerY))
  120. return (radius - distance(res)) / radius
  121. def _cv_small_disk(image_size):
  122. """Generates a disk level set function.
  123. The disk covers half of the image along its smallest dimension.
  124. """
  125. res = np.ones(image_size)
  126. centerY = int((image_size[0] - 1) / 2)
  127. centerX = int((image_size[1] - 1) / 2)
  128. res[centerY, centerX] = 0.0
  129. radius = float(min(centerX, centerY)) / 2.0
  130. return (radius - distance(res)) / (radius * 3)
  131. def _cv_init_level_set(init_level_set, image_shape, dtype=np.float64):
  132. """Generates an initial level set function conditional on input arguments."""
  133. if isinstance(init_level_set, str):
  134. if init_level_set == 'checkerboard':
  135. res = _cv_checkerboard(image_shape, 5, dtype)
  136. elif init_level_set == 'disk':
  137. res = _cv_large_disk(image_shape)
  138. elif init_level_set == 'small disk':
  139. res = _cv_small_disk(image_shape)
  140. else:
  141. raise ValueError("Incorrect name for starting level set preset.")
  142. else:
  143. res = init_level_set
  144. return res.astype(dtype, copy=False)
  145. def chan_vese(
  146. image,
  147. mu=0.25,
  148. lambda1=1.0,
  149. lambda2=1.0,
  150. tol=1e-3,
  151. max_num_iter=500,
  152. dt=0.5,
  153. init_level_set='checkerboard',
  154. extended_output=False,
  155. ):
  156. """Chan-Vese segmentation algorithm.
  157. Active contour model by evolving a level set. Can be used to
  158. segment objects without clearly defined boundaries.
  159. Parameters
  160. ----------
  161. image : (M, N) ndarray
  162. Grayscale image to be segmented.
  163. mu : float, optional
  164. 'edge length' weight parameter. Higher `mu` values will
  165. produce a 'round' edge, while values closer to zero will
  166. detect smaller objects.
  167. lambda1 : float, optional
  168. 'difference from average' weight parameter for the output
  169. region with value 'True'. If it is lower than `lambda2`, this
  170. region will have a larger range of values than the other.
  171. lambda2 : float, optional
  172. 'difference from average' weight parameter for the output
  173. region with value 'False'. If it is lower than `lambda1`, this
  174. region will have a larger range of values than the other.
  175. tol : float, positive, optional
  176. Level set variation tolerance between iterations. If the
  177. L2 norm difference between the level sets of successive
  178. iterations normalized by the area of the image is below this
  179. value, the algorithm will assume that the solution was
  180. reached.
  181. max_num_iter : uint, optional
  182. Maximum number of iterations allowed before the algorithm
  183. interrupts itself.
  184. dt : float, optional
  185. A multiplication factor applied at calculations for each step,
  186. serves to accelerate the algorithm. While higher values may
  187. speed up the algorithm, they may also lead to convergence
  188. problems.
  189. init_level_set : str or (M, N) ndarray, optional
  190. Defines the starting level set used by the algorithm.
  191. If a string is inputted, a level set that matches the image
  192. size will automatically be generated. Alternatively, it is
  193. possible to define a custom level set, which should be an
  194. array of float values, with the same shape as 'image'.
  195. Accepted string values are as follows.
  196. 'checkerboard'
  197. the starting level set is defined as
  198. sin(x/5*pi)*sin(y/5*pi), where x and y are pixel
  199. coordinates. This level set has fast convergence, but may
  200. fail to detect implicit edges.
  201. 'disk'
  202. the starting level set is defined as the opposite
  203. of the distance from the center of the image minus half of
  204. the minimum value between image width and image height.
  205. This is somewhat slower, but is more likely to properly
  206. detect implicit edges.
  207. 'small disk'
  208. the starting level set is defined as the
  209. opposite of the distance from the center of the image
  210. minus a quarter of the minimum value between image width
  211. and image height.
  212. extended_output : bool, optional
  213. If set to True, the return value will be a tuple containing
  214. the three return values (see below). If set to False which
  215. is the default value, only the 'segmentation' array will be
  216. returned.
  217. Returns
  218. -------
  219. segmentation : (M, N) ndarray, bool
  220. Segmentation produced by the algorithm.
  221. phi : (M, N) ndarray of floats
  222. Final level set computed by the algorithm.
  223. energies : list of floats
  224. Shows the evolution of the 'energy' for each step of the
  225. algorithm. This should allow to check whether the algorithm
  226. converged.
  227. Notes
  228. -----
  229. The Chan-Vese Algorithm is designed to segment objects without
  230. clearly defined boundaries. This algorithm is based on level sets
  231. that are evolved iteratively to minimize an energy, which is
  232. defined by weighted values corresponding to the sum of differences
  233. intensity from the average value outside the segmented region, the
  234. sum of differences from the average value inside the segmented
  235. region, and a term which is dependent on the length of the
  236. boundary of the segmented region.
  237. This algorithm was first proposed by Tony Chan and Luminita Vese,
  238. in a publication entitled "An Active Contour Model Without Edges"
  239. [1]_.
  240. This implementation of the algorithm is somewhat simplified in the
  241. sense that the area factor 'nu' described in the original paper is
  242. not implemented, and is only suitable for grayscale images.
  243. Typical values for `lambda1` and `lambda2` are 1. If the
  244. 'background' is very different from the segmented object in terms
  245. of distribution (for example, a uniform black image with figures
  246. of varying intensity), then these values should be different from
  247. each other.
  248. Typical values for mu are between 0 and 1, though higher values
  249. can be used when dealing with shapes with very ill-defined
  250. contours.
  251. The 'energy' which this algorithm tries to minimize is defined
  252. as the sum of the differences from the average within the region
  253. squared and weighed by the 'lambda' factors to which is added the
  254. length of the contour multiplied by the 'mu' factor.
  255. Supports 2D grayscale images only, and does not implement the area
  256. term described in the original article.
  257. References
  258. ----------
  259. .. [1] An Active Contour Model without Edges, Tony Chan and
  260. Luminita Vese, Scale-Space Theories in Computer Vision,
  261. 1999, :DOI:`10.1007/3-540-48236-9_13`
  262. .. [2] Chan-Vese Segmentation, Pascal Getreuer Image Processing On
  263. Line, 2 (2012), pp. 214-224,
  264. :DOI:`10.5201/ipol.2012.g-cv`
  265. .. [3] The Chan-Vese Algorithm - Project Report, Rami Cohen, 2011
  266. :arXiv:`1107.2782`
  267. """
  268. if len(image.shape) != 2:
  269. raise ValueError("Input image should be a 2D array.")
  270. float_dtype = _supported_float_type(image.dtype)
  271. phi = _cv_init_level_set(init_level_set, image.shape, dtype=float_dtype)
  272. if type(phi) != np.ndarray or phi.shape != image.shape:
  273. raise ValueError(
  274. "The dimensions of initial level set do not "
  275. "match the dimensions of image."
  276. )
  277. image = image.astype(float_dtype, copy=False)
  278. image = image - np.min(image)
  279. if np.max(image) != 0:
  280. image = image / np.max(image)
  281. i = 0
  282. old_energy = _cv_energy(image, phi, mu, lambda1, lambda2)
  283. energies = []
  284. phivar = tol + 1
  285. segmentation = phi > 0
  286. while phivar > tol and i < max_num_iter:
  287. # Save old level set values
  288. oldphi = phi
  289. # Calculate new level set
  290. phi = _cv_calculate_variation(image, phi, mu, lambda1, lambda2, dt)
  291. phi = _cv_reset_level_set(phi)
  292. phivar = np.sqrt(((phi - oldphi) ** 2).mean())
  293. # Extract energy and compare to previous level set and
  294. # segmentation to see if continuing is necessary
  295. segmentation = phi > 0
  296. new_energy = _cv_energy(image, phi, mu, lambda1, lambda2)
  297. # Save old energy values
  298. energies.append(old_energy)
  299. old_energy = new_energy
  300. i += 1
  301. if extended_output:
  302. return (segmentation, phi, energies)
  303. else:
  304. return segmentation