wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-07-07 08:48:01 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'neil_yager-greycomatrix'

Browse Source
This commit is contained in:
Stefan van der Walt
2011-11-30 09:49:32 -08:00
parent 4fbdeb4b41 917d1a80e8
commit 2db1b8367c
8 changed files with 565 additions and 3 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CONTRIBUTORS.txt
+1 -1
View File
@@ -79,7 +79,7 @@
Windows packaging and Python 3 compatibility.
- Neil Yager
Skeletonization.
Skeletonization and grey level co-occurrence matrices.
- Nelle Varoquaux
Renaming of the package to ``skimage``.
TASKS.txt
+1 -1
View File
@@ -47,7 +47,7 @@ modified to work as part of the scikit, others may be lacking in documentation
or tests.
* :strike:`Connected components`
* `Grey-level co-occurrence matrices <http://mentat.za.net/hg>`_
* :strike:`Grey-level co-occurrence matrices`
* Marching squares
* Nadav's bilateral filtering (first compare against CellProfiler's
code, based on http://groups.csail.mit.edu/graphics/bilagrid/bilagrid_web.pdf)
doc/examples/plot_glcm.py
+95
View File
@@ -0,0 +1,95 @@
"""
=====================
GLCM Texture Features
=====================
This example illustrates texture classification using texture
classification using grey level co-occurrence matrices (GLCMs).
A GLCM is a histogram of co-occurring greyscale values at a given
offset over an image.
In this example, samples of two different textures are extracted from
an image: grassy areas and sky areas. For each patch, a GLCM with
a horizontal offset of 5 is computed. Next, two features of the
GLCM matrices are computed: dissimilarity and correlation. These are
plotted to illustrate that the classes form clusters in feature space.
In a typical classification problem, the final step (not included in
this example) would be to train a classifier, such as logistic
regression, to label image patches from new images.
"""
from skimage.feature import greycomatrix, greycoprops
from skimage import data
import matplotlib.pyplot as plt
PATCH_SIZE = 21
# open the camera image
image = data.camera()
# select some patches from grassy areas of the image
grass_locations = [(474, 291), (440, 433), (466, 18), (462, 236)]
grass_patches = []
for loc in grass_locations:
grass_patches.append(image[loc[0]:loc[0] + PATCH_SIZE,
loc[1]:loc[1] + PATCH_SIZE])
# select some patches from sky areas of the image
sky_locations = [(54, 48), (21, 233), (90, 380), (195, 330)]
sky_patches = []
for loc in sky_locations:
sky_patches.append(image[loc[0]:loc[0] + PATCH_SIZE,
loc[1]:loc[1] + PATCH_SIZE])
# compute some GLCM properties each patch
xs = []
ys = []
for i, patch in enumerate(grass_patches + sky_patches):
glcm = greycomatrix(patch, [5], [0], 256, symmetric=True, normed=True)
xs.append(greycoprops(glcm, 'dissimilarity')[0, 0])
ys.append(greycoprops(glcm, 'correlation')[0, 0])
# create the figure
plt.figure(figsize=(8, 8))
# display the image patches
for i, patch in enumerate(grass_patches):
plt.subplot(3, len(grass_patches), len(grass_patches) * 1 + i + 1)
plt.imshow(patch, cmap=plt.cm.gray, interpolation='nearest',
vmin=0, vmax=255)
plt.xlabel('Grass %d' % (i + 1))
for i, patch in enumerate(sky_patches):
plt.subplot(3, len(grass_patches), len(grass_patches) * 2 + i + 1)
plt.imshow(patch, cmap=plt.cm.gray, interpolation='nearest',
vmin=0, vmax=255)
plt.xlabel('Sky %d' % (i + 1))
# display original image with locations of patches
plt.subplot(3, 2, 1)
plt.imshow(image, cmap=plt.cm.gray, interpolation='nearest',
vmin=0, vmax=255)
for (y, x) in grass_locations:
plt.plot(x + PATCH_SIZE / 2, y + PATCH_SIZE / 2, 'gs')
for (y, x) in sky_locations:
plt.plot(x + PATCH_SIZE / 2, y + PATCH_SIZE / 2, 'bs')
plt.xlabel('Original Image')
plt.xticks([])
plt.yticks([])
plt.axis('image')
# for each patch, plot (dissimilarity, correlation)
plt.subplot(3, 2, 2)
plt.plot(xs[:len(grass_patches)], ys[:len(grass_patches)], 'go',
label='Grass')
plt.plot(xs[len(grass_patches):], ys[len(grass_patches):], 'bo',
label='Sky')
plt.xlabel('GLCM Dissimilarity')
plt.ylabel('GLVM Correlation')
plt.legend()
# display the patches and plot
plt.suptitle('Grey level co-occurrence matrix features', fontsize=14)
plt.show()
skimage/feature/__init__.py
+2 -1
View File
@@ -1 +1,2 @@
from hog import hog
from hog import hog
from greycomatrix import greycomatrix, greycoprops
skimage/feature/_greycomatrix.pyx
+67
View File
@@ -0,0 +1,67 @@
"""Cython implementation for computing a grey level co-occurance matrix
"""
import numpy as np
cimport numpy as np
cimport cython
cdef extern from "math.h":
double sin(double)
double cos(double)
@cython.boundscheck(False)
def _glcm_loop(np.ndarray[dtype=np.uint8_t, ndim=2,
negative_indices=False, mode='c'] image,
np.ndarray[dtype=np.float64_t, ndim=1,
negative_indices=False, mode='c'] distances,
np.ndarray[dtype=np.float64_t, ndim=1,
negative_indices=False, mode='c'] angles,
int levels,
np.ndarray[dtype=np.uint32_t, ndim=4,
negative_indices=False, mode='c'] out
):
"""Perform co-occurnace matrix accumulation
Parameters
----------
image : ndarray
Input image, which is converted to the uint8 data type.
distances : ndarray
List of pixel pair distance offsets.
angles : ndarray
List of pixel pair angles in radians.
levels : int
The input image should contain integers in [0, levels-1],
where levels indicate the number of grey-levels counted
(typically 256 for an 8-bit image)
out : ndarray
On input a 4D array of zeros, and on output it contains
the results of the GLCM computation.
"""
cdef:
np.int32_t a_inx, d_idx
np.int32_t r, c, rows, cols, row, col
np.int32_t i, j
rows = image.shape[0]
cols = image.shape[1]
for a_idx, angle in enumerate(angles):
for d_idx, distance in enumerate(distances):
for r in range(rows):
for c in range(cols):
i = image[r, c]
# compute the location of the offset pixel
row = r + <int>(sin(angle) * distance + 0.5)
col = c + <int>(cos(angle) * distance + 0.5);
# make sure the offset is within bounds
if row >= 0 and row < rows and \
col >= 0 and col < cols:
j = image[row, col]
if i >= 0 and i < levels and \
j >= 0 and j < levels:
out[i, j, d_idx, a_idx] += 1
skimage/feature/greycomatrix.py
+225
View File
@@ -0,0 +1,225 @@
"""
Compute grey level co-occurrence matrices (GLCMs) and associated
properties to characterize image textures.
"""
import numpy as np
import skimage.util
from _greycomatrix import _glcm_loop
def greycomatrix(image, distances, angles, levels=256, symmetric=False,
normed=False):
"""Calculate the grey-level co-occurrence matrix.
A grey level co-occurence matrix is a histogram of co-occuring
greyscale values at a given offset over an image.
Parameters
----------
image : array_like of uint8
Integer typed input image. The image will be cast to uint8, so
the maximum value must be less than 256.
distances : array_like
List of pixel pair distance offsets.
angles : array_like
List of pixel pair angles in radians.
levels : int, optional
The input image should contain integers in [0, levels-1],
where levels indicate the number of grey-levels counted
(typically 256 for an 8-bit image). The maximum value is
256.
symmetric : bool, optional
If True, the output matrix `P[:, :, d, theta]` is symmetric. This
is accomplished by ignoring the order of value pairs, so both
(i, j) and (j, i) are accumulated when (i, j) is encountered
for a given offset. The default is False.
normed : bool, optional
If True, normalize each matrix `P[:, :, d, theta]` by dividing
by the total number of accumulated co-occurrences for the given
offset. The elements of the resulting matrix sum to 1. The
default is False.
Returns
-------
P : 4-D ndarray
The grey-level co-occurrence histogram. The value
`P[i,j,d,theta]` is the number of times that grey-level `j`
occurs at a distance `d` and at an angle `theta` from
grey-level `i`. If `normed` is `False`, the output is of
type uint32, otherwise it is float64.
References
----------
.. [1] The GLCM Tutorial Home Page,
http://www.fp.ucalgary.ca/mhallbey/tutorial.htm
.. [2] Pattern Recognition Engineering, Morton Nadler & Eric P.
Smith
.. [3] Wikipedia, http://en.wikipedia.org/wiki/Co-occurrence_matrix
Examples
--------
Compute 2 GLCMs: One for a 1-pixel offset to the right, and one
for a 1-pixel offset upwards.
>>> image = np.array([[0, 0, 1, 1],
... [0, 0, 1, 1],
... [0, 2, 2, 2],
... [2, 2, 3, 3]], dtype=np.uint8)
>>> result = greycomatrix(image, [1], [0, np.pi/2], levels=4)
>>> result[:, :, 0, 0]
array([[2, 2, 1, 0],
[0, 2, 0, 0],
[0, 0, 3, 1],
[0, 0, 0, 1]], dtype=uint32)
>>> result[:, :, 0, 1]
array([[3, 0, 2, 0],
[0, 2, 2, 0],
[0, 0, 1, 2],
[0, 0, 0, 0]], dtype=uint32)
"""
assert levels <= 256
image = np.ascontiguousarray(image)
assert image.ndim == 2
assert image.min() >= 0
assert image.max() < levels
image = image.astype(np.uint8)
distances = np.ascontiguousarray(distances, dtype=np.float64)
angles = np.ascontiguousarray(angles, dtype=np.float64)
assert distances.ndim == 1
assert angles.ndim == 1
P = np.zeros((levels, levels, len(distances), len(angles)),
dtype=np.uint32, order='C')
# count co-occurences
_glcm_loop(image, distances, angles, levels, P)
# make each GLMC symmetric
if symmetric:
P += np.transpose(P, (1, 0, 2, 3))
# normalize each GLMC
if normed:
P = P.astype(np.float64)
glcm_sums = np.apply_over_axes(np.sum, P, axes=(0, 1))
glcm_sums[glcm_sums == 0] = 1
P /= glcm_sums
return P
def greycoprops(P, prop='contrast'):
"""Calculate texture properties of a GLCM.
Compute a feature of a grey level co-occurrence matrix to serve as
a compact summary of the matrix. The properties are computed as
follows:
- 'contrast': :math:`\\sum_{i,j=0}^{levels-1} P_{i,j}(i-j)^2`
- 'dissimilarity': :math:`\\sum_{i,j=0}^{levels-1}P_{i,j}|i-j|`
- 'homogeneity': :math:`\\sum_{i,j=0}^{levels-1}\\frac{P_{i,j}}{1+(i-j)^2}`
- 'ASM': :math:`\\sum_{i,j=0}^{levels-1} P_{i,j}^2`
- 'energy': :math:`\\sqrt{ASM}`
- 'correlation':
.. math:: \\sum_{i,j=0}^{levels-1} P_{i,j}\\left[\\frac{(i-\\mu_i) \\
(j-\\mu_j)}{\\sqrt{(\\sigma_i^2)(\\sigma_j^2)}}\\right]
Parameters
----------
P : ndarray
Input array. `P` is the grey-level co-occurrence histogram
for which to compute the specified property. The value
`P[i,j,d,theta]` is the number of times that grey-level j
occurs at a distance d and at an angle theta from
grey-level i.
prop : {'contrast', 'dissimilarity', 'homogeneity', 'energy', \
'correlation', 'ASM'}, optional
The property of the GLCM to compute. The default is 'contrast'.
Returns
-------
results : 2-D ndarray
2-dimensional array. `results[d, a]` is the property 'prop' for
the d'th distance and the a'th angle.
References
----------
.. [1] The GLCM Tutorial Home Page,
http://www.fp.ucalgary.ca/mhallbey/tutorial.htm
Examples
--------
Compute the contrast for GLCMs with distances [1, 2] and angles
[0 degrees, 90 degrees]
>>> image = np.array([[0, 0, 1, 1],
... [0, 0, 1, 1],
... [0, 2, 2, 2],
... [2, 2, 3, 3]], dtype=np.uint8)
>>> g = greycomatrix(image, [1, 2], [0, np.pi/2], levels=4,
... normed=True, symmetric=True)
>>> contrast = greycoprops(g, 'contrast')
>>> contrast
array([[ 0.58333333, 1. ],
[ 1.25 , 2.75 ]])
"""
assert P.ndim == 4
(num_level, num_level2, num_dist, num_angle) = P.shape
assert num_level == num_level2
assert num_dist > 0
assert num_angle > 0
# create weights for specified property
I, J = np.ogrid[0:num_level, 0:num_level]
if prop == 'contrast':
weights = (I - J) ** 2
elif prop == 'dissimilarity':
weights = np.abs(I - J)
elif prop == 'homogeneity':
weights = 1. / (1. + (I - J) ** 2)
elif prop in ['ASM', 'energy', 'correlation']:
pass
else:
raise ValueError('%s is an invalid property' % (prop))
# compute property for each GLCM
if prop == 'energy':
asm = np.apply_over_axes(np.sum, (P ** 2), axes=(0, 1))[0, 0]
results = np.sqrt(asm)
elif prop == 'ASM':
results = np.apply_over_axes(np.sum, (P ** 2), axes=(0, 1))[0, 0]
elif prop == 'correlation':
results = np.zeros((num_dist, num_angle), dtype=np.float64)
I = np.array(range(num_level)).reshape((num_level, 1, 1, 1))
J = np.array(range(num_level)).reshape((1, num_level, 1, 1))
diff_i = I - np.apply_over_axes(np.sum, (I * P), axes=(0, 1))[0, 0]
diff_j = J - np.apply_over_axes(np.sum, (J * P), axes=(0, 1))[0, 0]
std_i = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_i) ** 2),
axes=(0, 1))[0, 0])
std_j = np.sqrt(np.apply_over_axes(np.sum, (P * (diff_j) ** 2),
axes=(0, 1))[0, 0])
cov = np.apply_over_axes(np.sum, (P * (diff_i * diff_j)),
axes=(0, 1))[0, 0]
# handle the special case of standard deviations near zero
mask_0 = std_i < 1e-15
mask_0[std_j < 1e-15] = True
results[mask_0] = 1
# handle the standard case
mask_1 = mask_0 == False
results[mask_1] = cov[mask_1] / (std_i[mask_1] * std_j[mask_1])
elif prop in ['contrast', 'dissimilarity', 'homogeneity']:
weights = weights.reshape((num_level, num_level, 1, 1))
results = np.apply_over_axes(np.sum, (P * weights), axes=(0, 1))[0, 0]
return results
skimage/feature/setup.py
+30
View File
@@ -0,0 +1,30 @@
#!/usr/bin/env python
import os
from skimage._build import cython
base_path = os.path.abspath(os.path.dirname(__file__))
def configuration(parent_package='', top_path=None):
from numpy.distutils.misc_util import Configuration, get_numpy_include_dirs
config = Configuration('feature', parent_package, top_path)
config.add_data_dir('tests')
cython(['_greycomatrix.pyx'], working_path=base_path)
config.add_extension('_greycomatrix', sources=['_greycomatrix.c'],
include_dirs=[get_numpy_include_dirs()])
return config
if __name__ == '__main__':
from numpy.distutils.core import setup
setup(maintainer = 'Scikits-image Developers',
author = 'Scikits-image Developers',
maintainer_email = 'scikits-image@googlegroups.com',
description = 'Features',
url = 'https://github.com/scikits-image/scikits-image',
license = 'SciPy License (BSD Style)',
**(configuration(top_path='').todict())
)
skimage/feature/tests/test_glcm.py
+144
View File
@@ -0,0 +1,144 @@
import numpy as np
from skimage.feature import greycomatrix, greycoprops
class TestGLCM():
def setup(self):
self.image = np.array([[0, 0, 1, 1],
[0, 0, 1, 1],
[0, 2, 2, 2],
[2, 2, 3, 3]], dtype=np.uint8)
def test_output_angles(self):
result = greycomatrix(self.image, [1], [0, np.pi / 2], 4)
assert result.shape == (4, 4, 1, 2)
expected1 = np.array([[2, 2, 1, 0],
[0, 2, 0, 0],
[0, 0, 3, 1],
[0, 0, 0, 1]], dtype=np.uint32)
np.testing.assert_array_equal(result[:, :, 0, 0], expected1)
expected2 = np.array([[3, 0, 2, 0],
[0, 2, 2, 0],
[0, 0, 1, 2],
[0, 0, 0, 0]], dtype=np.uint32)
np.testing.assert_array_equal(result[:, :, 0, 1], expected2)
def test_output_symmetric_1(self):
result = greycomatrix(self.image, [1], [np.pi / 2], 4,
symmetric=True)
assert result.shape == (4, 4, 1, 1)
expected = np.array([[6, 0, 2, 0],
[0, 4, 2, 0],
[2, 2, 2, 2],
[0, 0, 2, 0]], dtype=np.uint32)
np.testing.assert_array_equal(result[:, :, 0, 0], expected)
def test_output_distance(self):
im = np.array([[0, 0, 0, 0],
[1, 0, 0, 1],
[2, 0, 0, 2],
[3, 0, 0, 3]], dtype=np.uint8)
result = greycomatrix(im, [3], [0], 4, symmetric=False)
expected = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1]], dtype=np.uint32)
np.testing.assert_array_equal(result[:, :, 0, 0], expected)
def test_output_combo(self):
im = np.array([[0],
[1],
[2],
[3]], dtype=np.uint8)
result = greycomatrix(im, [1, 2], [0, np.pi / 2], 4)
assert result.shape == (4, 4, 2, 2)
z = np.zeros((4, 4), dtype=np.uint32)
e1 = np.array([[0, 1, 0, 0],
[0, 0, 1, 0],
[0, 0, 0, 1],
[0, 0, 0, 0]], dtype=np.uint32)
e2 = np.array([[0, 0, 1, 0],
[0, 0, 0, 1],
[0, 0, 0, 0],
[0, 0, 0, 0]], dtype=np.uint32)
np.testing.assert_array_equal(result[:, :, 0, 0], z)
np.testing.assert_array_equal(result[:, :, 1, 0], z)
np.testing.assert_array_equal(result[:, :, 0, 1], e1)
np.testing.assert_array_equal(result[:, :, 1, 1], e2)
def test_output_empty(self):
result = greycomatrix(self.image, [10], [0], 4)
np.testing.assert_array_equal(result[:, :, 0, 0],
np.zeros((4, 4), dtype=np.uint32))
result = greycomatrix(self.image, [10], [0], 4, normed=True)
np.testing.assert_array_equal(result[:, :, 0, 0],
np.zeros((4, 4), dtype=np.uint32))
def test_normed_symmetric(self):
result = greycomatrix(self.image, [1, 2, 3],
[0, np.pi / 2, np.pi], 4,
normed=True, symmetric=True)
for d in range(result.shape[2]):
for a in range(result.shape[3]):
np.testing.assert_almost_equal(result[:, :, d, a].sum(),
1.0)
np.testing.assert_array_equal(result[:, :, d, a],
result[:, :, d, a].transpose())
def test_contrast(self):
result = greycomatrix(self.image, [1, 2], [0], 4,
normed=True, symmetric=True)
result = np.round(result, 3)
contrast = greycoprops(result, 'contrast')
np.testing.assert_almost_equal(contrast[0, 0], 0.586)
def test_dissimilarity(self):
result = greycomatrix(self.image, [1], [0, np.pi / 2], 4,
normed=True, symmetric=True)
result = np.round(result, 3)
dissimilarity = greycoprops(result, 'dissimilarity')
np.testing.assert_almost_equal(dissimilarity[0, 0], 0.418)
def test_dissimilarity_2(self):
result = greycomatrix(self.image, [1, 3], [np.pi/2], 4,
normed=True, symmetric=True)
result = np.round(result, 3)
dissimilarity = greycoprops(result, 'dissimilarity')[0, 0]
np.testing.assert_almost_equal(dissimilarity, 0.664)
def test_invalid_property(self):
result = greycomatrix(self.image, [1], [0], 4)
np.testing.assert_raises(ValueError, greycoprops,
result, 'ABC')
def test_homogeneity(self):
result = greycomatrix(self.image, [1], [0, 6], 4, normed=True,
symmetric=True)
homogeneity = greycoprops(result, 'homogeneity')[0, 0]
np.testing.assert_almost_equal(homogeneity, 0.80833333)
def test_energy(self):
result = greycomatrix(self.image, [1], [0, 4], 4, normed=True,
symmetric=True)
energy = greycoprops(result, 'energy')[0, 0]
np.testing.assert_almost_equal(energy, 0.38188131)
def test_correlation(self):
result = greycomatrix(self.image, [1, 2], [0], 4, normed=True,
symmetric=True)
energy = greycoprops(result, 'correlation')
np.testing.assert_almost_equal(energy[0, 0], 0.71953255)
np.testing.assert_almost_equal(energy[1, 0], 0.41176470)
def test_uniform_properties(self):
im = np.ones((4, 4), dtype=np.uint8)
result = greycomatrix(im, [1, 2, 8], [0, np.pi / 2], 4, normed=True,
symmetric=True)
for prop in ['contrast', 'dissimilarity', 'homogeneity',
'energy', 'correlation', 'ASM']:
greycoprops(result, prop)
if __name__ == '__main__':
np.testing.run_module_suite()