quantimpy.segmentation module¶
Functions for image segmentation
This module contains various functions for the image segmentation of both 2D and 3D NumPy 1 arrays.
- quantimpy.segmentation.anisodiff(image, option, niter, K, gamma)¶
Anisotropic diffusion filter
This function applies an anisotropic diffusion filter to the 2D or 3D NumPy array image. This is also known as Perona Malik diffusion 2 and is an edge preserving noise reduction method. The code is based on a Matlab code by Perona, Shiota, and Malik 3.
- Parameters
image (ndarray, {int, uint, float}) – 2D or 3D grayscale input image.
option (int, defaults to 1) –
The option parameter selects the conduction coefficient used by the filter. option=0 selects the following conduction coefficient:
\[g (\nabla I) = \exp{(-\frac{||\nabla I||}{K})},\]where \(\nabla I\) is the image brightness gradient, and \(K\) is a constant. This equation is used in a Matlab code by Perona, Shiota, and Malik 3. option=1 selects the conduction coefficient:
\[g (\nabla I) = \exp{(-\left(\frac{||\nabla I||}{K}\right)^{2})},\]and option=2 selects the coefficient:
\[g (\nabla I) = \frac{1}{1 + (\frac{||\nabla I||}{K})^{2}}.\]Option one privileges high-contrast edges over low-contrast ones, while option two privileges wide regions over smaller ones 2.
niter (int, defaults to 1) – The number of iterations that the filter is applied.
K (float, defaults to 50) – The value of constant \(K\) in the above equations.
gamma (float, defaults to 0.1) – Sets the diffusion “time” step size. When \(\gamma \leq 0.25\), stability is ensured.
- Returns
out – Noise reduced 2D or 3D output image. The return data type is float and the image is normalized betweeen 0 and 1 or -1 and 1.
- Return type
ndarray, float
Examples
This example uses the NumPy 1, and Matplotlib Python packages 4. The NumPy data file “rock_2d.npy” is available on Github 9 10.
import numpy as np import matplotlib.pyplot as plt from quantimpy import segmentation as sg # Load data image = np.load("rock_2d.npy") # Apply anisotropic diffusion filter diffusion = sg.anisodiff(image, niter=5) # Show results fig = plt.figure() plt.gray() ax1 = fig.add_subplot(121) # left side ax2 = fig.add_subplot(122) # right side ax1.imshow(image) ax2.imshow(diffusion) plt.show()
- quantimpy.segmentation.bilevel(image, thres_min, thres_max, debug=None)¶
Compute bi-level image segmentation
This function computes the bi-level binary segmentation of the 2D or 3D NumPy array image 7. First, an initial segmentation is computed using the minimum threshold value thres_min. Then, iteratively, this initial segmented image is dilated and all voxels that are below the maximum threshold value thres_max are added to the segmented image. This continues untill no more changes to the segmented image are observed.
- Parameters
image (ndarray, {int, uint, float}) – 2D or 3D grayscale input image.
thres_min (float) – Minimum threshold value for bi-level segmentation.
thres_max (float) – Maximum threshold value for bi-level segmentation.
debug (str, defaults to "None") – Output directory for debugging images. When this parameter is set, filtered images, masks, and histrograms are written to disk. Set to “./” to write to the working directory. The default is “None”.
- Returns
out – Returns either a 2D or 3D boolean ndarray of the segmented image.
- Return type
ndarray, bool
See also
Examples
This example uses the NumPy 1, and Matplotlib Python packages 4. The NumPy data file “rock_3d.npy” is available on Github 9 10.
import numpy as np import matplotlib.pyplot as plt from quantimpy import segmentation as sg # Load data image = np.load("rock_3d.npy") # Apply anisotropic diffusion filter diffusion = sg.anisodiff(image, niter=3) # Compute minimum and maximum thresholds thrshld_min, thrshld_max = sg.gradient(diffusion, alpha=1.5) # Apply bi-level segmentation binary = sg.bilevel(diffusion, thrshld_min, thrshld_max) # Show results fig = plt.figure() plt.gray() # show the result in grayscale ax1 = fig.add_subplot(131) ax2 = fig.add_subplot(132) ax3 = fig.add_subplot(133) ax1.imshow(image[50,:,:]) ax2.imshow(diffusion[50,:,:]) ax3.imshow(binary[50,:,:]) plt.show()
- quantimpy.segmentation.gradient(image, alpha=1.25, debug=None)¶
Compute thresholds for bi-level segmentation using gradient masks
This function computes the minimum and maximum threshold values for the 2D or 3D NumPy array image using gradient masks 8. These threshold values, in turn, can be used as inputs for the bi-level segmentation function,
bilevel(). The gradient masks are computed using Sobel and Laplace edge detection filters combined with the unimodal thresholding function,unimodal().- Parameters
image (ndarray, {int, uint, float}) – 2D or 3D grayscale input image.
alpha (float, defaults to 1.25) –
The parameter \(\alpha > 1\) is used to compute the minimum threshold value using the formula:
\[T_{\text{min}} = x_{\text{mode}} - \alpha (x_{\text{mode}} - T_{\text{max}}),\]where \(T_{\text{min}}\) is the minimum threshold value, \(x_{\text{mode}}\) is the mode of the histogram of image, and \(T_{\text{max}}\) is the maximum threshold value. The less noise image contains, the closer \(\alpha\) can be set to one.
debug (str, defaults to "None") – Output directory for debugging images. When this parameter is set, filtered images, masks, and histrograms are written to disk. Set to “./” to write to the working directory. The default is “None”.
- Returns
out – Returns minimum and maximum threshold values
- Return type
tuple, float
See also
Examples
This example uses the NumPy package 1. The NumPy data file “rock_2d.npy” is available on Github 9 10.
import numpy as np from quantimpy import segmentation as sg # Load data image = np.load("rock_2d.npy") # Apply anisotropic diffusion filter diffusion = sg.anisodiff(image, niter=5) # Compute minimum and maximum thresholds thrshld_min, thrshld_max = sg.gradient(diffusion) # Print results print(thrshld_min, thrshld_max)
- quantimpy.segmentation.histogram(image, bits=8)¶
Create an image histogram
This function creates an histogram for the 2D or 3D NumPy array image. The histogram is 8-bit (\(2^8\) bins) by default. The function is coded around the numpy.histogram function. However, this functions returns the center locations of the bins instead of the edges. For float or 16-bit images the bin size is scaled accordingly.
- Parameters
image (ndarray, {int, uint, float}) – 2D or 3D grayscale input image.
bits (int, defaults to 8) – \(2^{\text{bits}}\) bins are used for the histogram. Defaults to 8 bits or 256 bins.
- Returns
out – Returns two ndarrays: one with the histogram and one with the bin centers.
- Return type
tuple, float
See also
Examples
This example uses the NumPy 1, and Matplotlib Python packages 4. The NumPy data file “rock_3d.npy” is available on Github 9 10.
import numpy as np import matplotlib.pyplot as plt from quantimpy import segmentation as sg # Load data uint16 image = np.load("rock_3d.npy") # Compute histpgram hist, bins = sg.histogram(image, bits=8) width = bins[1] - bins[0] print(bins[0],bins[-1]) # Plot histogram plt.bar(bins,hist,width=width) plt.show()
- quantimpy.segmentation.unimodal(hist, side='right')¶
Compute unimodal threshold
Using image histogram hist, this function computes the unimodal threshold 6. This algorithms is slightly modified modified from the original method. Instead of defining the end of the distribution as the point where the histogram is zero, this algorithm takes the point that contains 99.7% of the observations. This is equivalent to three times the standard deviation for a Gaussian distribution.
- Parameters
hist (ndarray, int) – Histogram computed by the function
histogram().side (string, defaults to "right") – Whether to compute the unimodal threshold on the left or right side of the histogram maximum. Defaults to “right”.
- Returns
out – Index of the unimodal threshold value in input array hist.
- Return type
int
See also
Examples
This example uses the NumPy 1, Matplotlib 4, and SciPy Python packages 5. NumPy data file “rock_2d.npy” is available on Github 9 10.
import numpy as np import matplotlib.pyplot as plt from scipy import ndimage from quantimpy import filters # Load data image = np.load("rock_2d.npy") # Filter image result = filters.anisodiff(image, niter=3) # Edge detection laplace = ndimage.laplace(result) # Compute histpgram hist, bins = filters.histogram(laplace) width = bins[1] - bins[0] # Compute unimodal threshold thrshld = filters.unimodal(hist) # Plot histogram plt.bar(bins, hist,width=width) plt.scatter(bins[thrshld], hist[thrshld]) plt.show() # Compute unimodal threshold thrshld = filters.unimodal(hist, side="left") # Plot histogram plt.bar(bins, hist, width=width) plt.scatter(bins[thrshld], hist[thrshld]) plt.show()
References
- 1(1,2,3,4,5,6)
Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt et al., “Array programming with NumPy”, Nature, vol. 585, pp 357-362, 2020, doi:10.1038/s41586-020-2649-2
- 2(1,2)
Pietro Perona and Jitendra Malik, “Scale-space and edge detection using anisotropic diffusion”, IEEE Transactions on pattern analysis and machine intelligence, vol. 12, no. 7, pp 629-639, 1990, doi:10.1109/34.56205
- 3(1,2)
Pietro Perona, Takahiro Shiota, and Jitendra Malik, “Anisotropic diffusion”, in “Geometry-driven diffusion in computer vision”, ed. Bart M. ter Haar Romeny, pp 73-92, 1994, isbn: 9789401716994
- 4(1,2,3,4)
John D. Hunter, “Matplotlib: A 2D Graphics Environment”, Computing in Science & Engineering, vol. 9, no. 3, pp. 90-95, 2007. doi:10.1109/MCSE.2007.55
- 5
Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, et al., “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python”, Nature Methods, vol. 17, pp 261-272, 2020, doi:10.1038/s41592-019-0686-2
- 6
Paul Rosin, “Unimodal thresholding”, Pattern recognition, vol. 34, no. 11, pp 2083-2096, 2001, doi:10.1016/S0031-3203(00)00136-9
- 7
Hans-Jörg Vogel and Andre Kretzschmar, “Topological characterization of pore space in soil—sample preparation and digital image-processing”, Geoderma, vol. 73, no. 1-2, pp 23–38, 1996, doi:10.1016/0016-7061(96)00043-2
- 8
Steffen Schlüter, Ulrich Weller, and Hans-Jörg Vogel, “Segmentation of X-ray microtomography images of soil using gradient masks”, Computers & Geosciences, vol. 36, no. 10, pp 1246–1251, 2010, doi:10.1016/j.cageo.2010.02.007