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

Merge pull request #239 from ahojnnes/texture

Browse Source 
ENH: Add local binary patterns for texture characterization.

Conflicts:
	CONTRIBUTORS.txt
This commit is contained in:
Stefan van der Walt
2012-08-20 14:37:30 -07:00
parent 56cbdf6a01 014b490578
commit 73332f9959
14 changed files with 699 additions and 149 deletions
Download Patch File Download Diff File Expand all files Collapse all files
CONTRIBUTORS.txt
+2 -1
View File
@@ -108,6 +108,7 @@
Adaptive thresholding
Implementation of Matlab's `regionprops`
Estimation of geometric transformation parameters
Local binary pattern texture classification
- Pavel Campr
Fixes and tests for Histograms of Oriented Gradients.
doc/examples/plot_local_binary_pattern.py
+87
View File
@@ -0,0 +1,87 @@
"""
===============================================
Local Binary Pattern for texture classification
===============================================
In this example, we will see how to classify textures based on LBP (Local
Binary Pattern). The histogram of the LBP result is a good measure to classify
textures. For simplicity the histogram distributions are then tested against
each other using the Kullback-Leibler-Divergence.
"""
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import scipy.ndimage as nd
import skimage.feature as ft
from skimage import data
# settings for LBP
METHOD = 'uniform'
P = 16
R = 2
matplotlib.rcParams['font.size'] = 9
def kullback_leibler_divergence(p, q):
p = np.asarray(p)
q = np.asarray(q)
filt = np.logical_and(p != 0, q != 0)
return np.sum(p[filt] * np.log2(p[filt] / q[filt]))
def match(refs, img):
best_score = 10
best_name = None
lbp = ft.local_binary_pattern(img, P, R, METHOD)
hist, _ = np.histogram(lbp, normed=True, bins=P + 2, range=(0, P + 2))
for name, ref in refs.items():
ref_hist, _ = np.histogram(ref, normed=True, bins=P + 2,
range=(0, P + 2))
score = kullback_leibler_divergence(hist, ref_hist)
if score < best_score:
best_score = score
best_name = name
return best_name
brick = data.load('brick.png')
grass = data.load('grass.png')
wall = data.load('rough-wall.png')
refs = {
'brick': ft.local_binary_pattern(brick, P, R, METHOD),
'grass': ft.local_binary_pattern(grass, P, R, METHOD),
'wall': ft.local_binary_pattern(wall, P, R, METHOD)
}
# classify rotated textures
print 'Rotated images matched against references using LBP:'
print 'original: brick, rotated: 30deg, match result:',
print match(refs, nd.rotate(brick, angle=30, reshape=False))
print 'original: brick, rotated: 70deg, match result:',
print match(refs, nd.rotate(brick, angle=70, reshape=False))
print 'original: grass, rotated: 145deg, match result:',
print match(refs, nd.rotate(grass, angle=145, reshape=False))
# plot histograms of LBP of textures
fig, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(nrows=2, ncols=3,
figsize=(9, 6))
plt.gray()
ax1.imshow(brick)
ax1.axis('off')
ax4.hist(refs['brick'].ravel(), normed=True, bins=P + 2, range=(0, P + 2))
ax4.set_ylabel('Percentage')
ax2.imshow(grass)
ax2.axis('off')
ax5.hist(refs['grass'].ravel(), normed=True, bins=P + 2, range=(0, P + 2))
ax5.set_xlabel('Uniform LBP values')
ax3.imshow(wall)
ax3.axis('off')
ax6.hist(refs['wall'].ravel(), normed=True, bins=P + 2, range=(0, P + 2))
plt.show()
skimage/data/brick.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 170 KiB

skimage/data/grass.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 199 KiB

skimage/data/rough-wall.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 193 KiB

skimage/feature/__init__.py
+1 -1
View File
@@ -1,5 +1,5 @@
from ._hog import hog
from ._greycomatrix import greycomatrix, greycoprops
from .texture import greycomatrix, greycoprops, local_binary_pattern
from .peak import peak_local_max
from ._harris import harris
from .template import match_template
skimage/feature/_greycomatrix_cy.pyx
-67
View File
@@ -1,67 +0,0 @@
"""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/_texture.pyx
+181
View File
@@ -0,0 +1,181 @@
#cython: cdivison=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
import numpy as np
cimport numpy as np
from libc.math cimport sin, cos, abs
from skimage.transform._project cimport bilinear_interpolation
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-occurrence 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
cdef inline int _bit_rotate_right(int value, int length):
"""Cyclic bit shift to the right.
Parameters
----------
value : int
integer value to shift
length : int
number of bits of integer
"""
return (value >> 1) | ((value & 1) << (length - 1))
def _local_binary_pattern(np.ndarray[double, ndim=2] image,
int P, float R, char method='D'):
"""Gray scale and rotation invariant LBP (Local Binary Patterns).
LBP is an invariant descriptor that can be used for texture classification.
Parameters
----------
image : (N, M) double array
Graylevel image.
P : int
Number of circularly symmetric neighbour set points (quantization of the
angular space).
R : float
Radius of circle (spatial resolution of the operator).
method : {'D', 'R', 'U', 'V'}
Method to determine the pattern::
* 'D': 'default'
* 'R': 'ror'
* 'U': 'uniform'
* 'V': 'var'
Returns
-------
output : (N, M) array
LBP image.
"""
# texture weights
cdef np.ndarray[int, ndim=1] weights = 2 ** np.arange(P, dtype=np.int32)
# local position of texture elements
rp = - R * np.sin(2 * np.pi * np.arange(P, dtype=np.double) / P)
cp = R * np.cos(2 * np.pi * np.arange(P, dtype=np.double) / P)
cdef np.ndarray[double, ndim=2] coords = np.round(np.vstack([rp, cp]).T, 5)
# pre allocate arrays for computation
cdef np.ndarray[double, ndim=1] texture = np.zeros(P, np.double)
cdef np.ndarray[char, ndim=1] signed_texture = np.zeros(P, np.int8)
cdef np.ndarray[int, ndim=1] rotation_chain = np.zeros(P, np.int32)
output_shape = (image.shape[0], image.shape[1])
cdef np.ndarray[double, ndim=2] output = np.zeros(output_shape, np.double)
cdef int rows = image.shape[0]
cdef int cols = image.shape[1]
cdef double lbp
cdef int r, c, changes, i
for r in range(image.shape[0]):
for c in range(image.shape[1]):
for i in range(P):
texture[i] = bilinear_interpolation(<double*>image.data,
rows, cols, r + coords[i, 0], c + coords[i, 1], 'C')
# signed / thresholded texture
for i in range(P):
if texture[i] - image[r, c] >= 0:
signed_texture[i] = 1
else:
signed_texture[i] = 0
lbp = 0
# if method == 'uniform' or method == 'var':
if method == 'U' or method == 'V':
# determine number of 0 - 1 changes
changes = 0
for i in range(P - 1):
changes += abs(signed_texture[i] - signed_texture[i + 1])
if changes <= 2:
for i in range(P):
lbp += signed_texture[i]
else:
lbp = P + 1
if method == 'V':
var = np.var(texture)
if var != 0:
lbp /= var
else:
lbp = np.nan
else:
# method == 'default'
for i in range(P):
lbp += signed_texture[i] * weights[i]
# method == 'ror'
if method == 'R':
# shift LBP P times to the right and get minimum value
rotation_chain[0] = <int>lbp
for i in range(1, P):
rotation_chain[i] = \
_bit_rotate_right(rotation_chain[i - 1], P)
lbp = rotation_chain[0]
for i in range(1, P):
lbp = min(lbp, rotation_chain[i])
output[r, c] = lbp
return output
skimage/feature/setup.py
+4 -3
View File
@@ -12,11 +12,12 @@ def configuration(parent_package='', top_path=None):
config = Configuration('feature', parent_package, top_path)
config.add_data_dir('tests')
cython(['_greycomatrix_cy.pyx'], working_path=base_path)
cython(['_texture.pyx'], working_path=base_path)
cython(['_template.pyx'], working_path=base_path)
config.add_extension('_greycomatrix_cy', sources=['_greycomatrix_cy.c'],
include_dirs=[get_numpy_include_dirs()])
config.add_extension('_texture', sources=['_texture.c'],
include_dirs=[get_numpy_include_dirs(),
'../transform'])
config.add_extension('_template', sources=['_template.c'],
include_dirs=[get_numpy_include_dirs()])
skimage/feature/tests/test_glcm.py → skimage/feature/tests/test_texture.py
+60 -1
View File
@@ -1,8 +1,9 @@
import numpy as np
from skimage.feature import greycomatrix, greycoprops
from skimage.feature import greycomatrix, greycoprops, local_binary_pattern
class TestGLCM():
def setup(self):
self.image = np.array([[0, 0, 1, 1],
[0, 0, 1, 1],
@@ -140,5 +141,63 @@ class TestGLCM():
'energy', 'correlation', 'ASM']:
greycoprops(result, prop)
class TestLBP():
def setup(self):
self.image = np.array([[255, 6, 255, 0, 141, 0],
[ 48, 250, 204, 166, 223, 63],
[ 8, 0, 159, 50, 255, 30],
[167, 255, 63, 40, 128, 255],
[ 0, 255, 30, 34, 255, 24],
[146, 241, 255, 0, 189, 126]], dtype='double')
def test_default(self):
lbp = local_binary_pattern(self.image, 8, 1, 'default')
ref = np.array([[ 0, 251, 0, 255, 96, 255],
[143, 0, 20, 153, 64, 56],
[238, 255, 12, 191, 0, 252],
[129, 64., 62, 159, 199, 0],
[255, 4, 255, 175, 0, 254],
[ 3, 5, 0, 255, 4, 24]])
np.testing.assert_array_equal(lbp, ref)
def test_ror(self):
lbp = local_binary_pattern(self.image, 8, 1, 'ror')
ref = np.array([[ 0, 127, 0, 255, 3, 255],
[ 31, 0, 5, 51, 1, 7],
[119, 255, 3, 127, 0, 63],
[ 3, 1, 31, 63, 31, 0],
[255, 1, 255, 95, 0, 127],
[ 3, 5, 0, 255, 1, 3]])
np.testing.assert_array_equal(lbp, ref)
def test_uniform(self):
lbp = local_binary_pattern(self.image, 8, 1, 'uniform')
ref = np.array([[0, 7, 0, 8, 2, 8],
[5, 0, 9, 9, 1, 3],
[9, 8, 2, 7, 0, 6],
[2, 1, 5, 6, 5, 0],
[8, 1, 8, 9, 0, 7],
[2, 9, 0, 8, 1, 2]])
np.testing.assert_array_equal(lbp, ref)
def test_var(self):
lbp = local_binary_pattern(self.image, 8, 1, 'var')
ref = np.array([[0. , 0.00072786, 0. , 0.00115377,
0.00032355, 0.00224467],
[0.00051758, 0. , 0.0026383 , 0.00163246,
0.00027414, 0.00041124],
[0.00192834, 0.00130368, 0.00042095, 0.00171894,
0. , 0.00063726],
[0.00023048, 0.00019464 , 0.00082291, 0.00225386,
0.00076696, 0. ],
[0.00097253, 0.00013236, 0.0009134 , 0.0014467 ,
0. , 0.00082472],
[0.00024701, 0.0012277 , 0. , 0.00109869,
0.00015445, 0.00035881]])
np.testing.assert_array_almost_equal(lbp, ref)
if __name__ == '__main__':
np.testing.run_module_suite()
skimage/feature/texture.py
+278
View File
@@ -0,0 +1,278 @@
"""
Methods to characterize image textures.
"""
import math
import numpy as np
from scipy import ndimage
from ._texture import _glcm_loop, _local_binary_pattern
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:
Pt = np.transpose(P, (1, 0, 2, 3))
P = P + Pt
# 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
def local_binary_pattern(image, P, R, method='default'):
"""Gray scale and rotation invariant LBP (Local Binary Patterns).
LBP is an invariant descriptor that can be used for texture classification.
Parameters
----------
image : (N, M) array
Graylevel image.
P : int
Number of circularly symmetric neighbour set points (quantization of the
angular space).
R : float
Radius of circle (spatial resolution of the operator).
method : {'D', 'R', 'U', 'V'}
Method to determine the pattern::
* 'default': original local binary pattern which is gray scale but not
rotation invariant.
* 'ror': extension of default implementation which is gray scale and
rotation invariant.
* 'uniform': improved rotation invariance with uniform patterns and
finer quantization of the angular space which is gray scale and
rotation invariant.
* 'var': rotation invariant variance measures of the contrast of local
image texture which is rotation but not gray scale invariant.
Returns
-------
output : (N, M) array
LBP image.
References
----------
.. [1] Multiresolution Gray-Scale and Rotation Invariant Texture
Classification with Local Binary Patterns.
Timo Ojala, Matti Pietikainen, Topi Maenpaa.
http://www.rafbis.it/biplab15/images/stories/docenti/Danielriccio/\
Articoliriferimento/LBP.pdf, 2002.
"""
methods = {
'default': ord('D'),
'ror': ord('R'),
'uniform': ord('U'),
'var': ord('V')
}
image = np.array(image, dtype='double', copy=True)
output = _local_binary_pattern(image, P, R, methods[method.lower()])
return output
skimage/transform/_project.pxd
+7
View File
@@ -0,0 +1,7 @@
cimport numpy as np
import numpy as np
cdef inline double bilinear_interpolation(double* image, int rows, int cols,
double r, double c, char mode,
double cval=*)
skimage/transform/_project.pyx
+68 -64
View File
@@ -1,38 +1,66 @@
#cython: cdivison=True boundscheck=False
#cython: cdivison=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
__all__ = ['homography']
cimport cython
cimport numpy as np
import numpy as np
import cython
from cython.operator import dereference
from libc.math cimport ceil, floor
np.import_array()
cdef extern from "math.h":
double floor(double)
double fmod(double, double)
cdef inline double bilinear_interpolation(double* image, int rows, int cols,
double r, double c, char mode,
double cval=0):
"""Bilinear interpolation at a given position in the image.
cdef double get_pixel(double *image, int rows, int cols,
int r, int c, char mode, double cval=0):
Parameters
----------
image : double array
Input image.
rows, cols: int
Shape of image.
r, c : int
Position at which to interpolate.
mode : {'C', 'W', 'M'}
Wrapping mode. Constant, Wrap or Mirror.
cval : double
Constant value to use for constant mode.
"""
cdef double dr, dc
cdef int minr, minc, maxr, maxc
minr = <int>floor(r)
minc = <int>floor(c)
maxr = <int>ceil(r)
maxc = <int>ceil(c)
dr = r - minr
dc = c - minc
top = (1 - dc) * get_pixel(image, rows, cols, minr, minc, mode, cval) \
+ dc * get_pixel(image, rows, cols, minr, maxc, mode, cval)
bottom = (1 - dc) * get_pixel(image, rows, cols, maxr, minc, mode, cval) \
+ dc * get_pixel(image, rows, cols, maxr, maxc, mode, cval)
return (1 - dr) * top + dr * bottom
cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
char mode, double cval=0):
"""Get a pixel from the image, taking wrapping mode into consideration.
Parameters
----------
image : *double
image : double array
Input image.
rows, cols : int
Dimensions of image.
rows, cols: int
Shape of image.
r, c : int
Position at which to get the pixel.
mode : {'C', 'W', 'M'}
Wrapping mode. Constant, Wrap or Mirror.
Wrapping mode. Constant, Wrap or Mirror.
cval : double
Constant value to use for mode constant.
Constant value to use for constant mode.
"""
if mode == 'C':
if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1):
@@ -40,13 +68,13 @@ cdef double get_pixel(double *image, int rows, int cols,
else:
return image[r * cols + c]
else:
return image[coord_map(rows, r, mode) * cols +
coord_map(cols, c, mode)]
return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)]
cdef int coord_map(int dim, int coord, char mode):
cdef inline int coord_map(int dim, int coord, char mode):
"""
Wrap a coordinate, according to a given dimension and mode.
Wrap a coordinate, according to a given mode.
Parameters
----------
dim : int
@@ -56,7 +84,7 @@ cdef int coord_map(int dim, int coord, char mode):
mode : {'W', 'M'}
Whether to wrap or mirror the coordinate if it
falls outside [0, dim).
"""
dim = dim - 1
if mode == 'M': # mirror
@@ -79,7 +107,9 @@ cdef int coord_map(int dim, int coord, char mode):
return coord
cdef tf(double x, double y, double* H, double *x_, double *y_):
cdef inline _matrix_transform(double x, double y, double* H, double *x_,
double *y_):
"""Apply a homography to a coordinate.
Parameters
@@ -98,18 +128,15 @@ cdef tf(double x, double y, double* H, double *x_, double *y_):
yy = H[3] * x + H[4] * y + H[5]
zz = H[6] * x + H[7] * y + H[8]
xx = xx / zz
yy = yy / zz
x_[0] = xx / zz
y_[0] = yy / zz
x_[0] = xx
y_[0] = yy
@cython.boundscheck(False)
def homography(np.ndarray image, np.ndarray H, output_shape=None,
mode='constant', double cval=0):
"""
Projective transformation (homography).
Perform a projective transformation (homography) of a
floating point image, using bi-linear interpolation.
@@ -140,8 +167,6 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
Transformation matrix H that defines the homography.
output_shape : tuple (rows, cols)
Shape of the output image generated.
order : int
Order of splines used in interpolation.
mode : {'constant', 'mirror', 'wrap'}
How to handle values outside the image borders.
cval : string
@@ -150,8 +175,7 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
"""
cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] img = \
np.ascontiguousarray(image, dtype=np.double)
cdef np.ndarray[dtype=np.double_t, ndim=2] img = image.astype(np.double)
cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] M = \
np.ascontiguousarray(np.linalg.inv(H))
@@ -165,7 +189,6 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
elif mode == 'mirror':
mode_c = ord('M')
cdef int out_r, out_c, columns, rows
if output_shape is None:
out_r = img.shape[0]
out_c = img.shape[1]
@@ -173,37 +196,18 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
out_r = output_shape[0]
out_c = output_shape[1]
rows = img.shape[0]
columns = img.shape[1]
cdef np.ndarray[dtype=np.double_t, ndim=2] out = \
np.zeros((out_r, out_c), dtype=np.double)
cdef int tfr, tfc, r_int, c_int
cdef double y0, y1, y2, y3
cdef double r, c, z, t, u
cdef int tfr, tfc
cdef double r, c
cdef int rows = img.shape[0]
cdef int cols = img.shape[1]
for tfr in range(out_r):
for tfc in range(out_c):
tf(tfc, tfr, <double*>M.data, &c, &r)
r_int = <int>floor(r)
c_int = <int>floor(c)
t = r - r_int
u = c - c_int
y0 = get_pixel(<double*>img.data, rows, columns,
r_int, c_int, mode_c)
y1 = get_pixel(<double*>img.data, rows, columns,
r_int + 1, c_int, mode_c)
y2 = get_pixel(<double*>img.data, rows, columns,
r_int + 1, c_int + 1, mode_c)
y3 = get_pixel(<double*>img.data, rows, columns,
r_int, c_int + 1, mode_c)
out[tfr, tfc] = \
(1 - t) * (1 - u) * y0 + \
t * (1 - u) * y1 + \
t * u * y2 + (1 - t) * u * y3;
_matrix_transform(tfc, tfr, <double*>M.data, &c, &r)
out[tfr, tfc] = bilinear_interpolation(<double*>img.data, rows,
cols, r, c, mode_c)
return out
skimage/transform/tests/test_warps.py
+11 -12
View File
@@ -2,7 +2,7 @@ from numpy.testing import assert_array_almost_equal, run_module_suite
import numpy as np
from skimage.transform import (warp, homography, fast_homography,
SimilarityTransform)
SimilarityTransform, ProjectiveTransform)
from skimage import transform as tf, data, img_as_float
from skimage.color import rgb2gray
@@ -34,7 +34,7 @@ def test_homography():
def test_fast_homography():
img = rgb2gray(data.lena()).astype(np.uint8)
img = rgb2gray(data.lena())
img = img[:, :100]
theta = np.deg2rad(30)
@@ -49,20 +49,19 @@ def test_fast_homography():
H[:2, 2] = [tx, ty]
for mode in ('constant', 'mirror', 'wrap'):
p0 = homography(img, H, mode=mode, order=1)
p0 = warp(img, ProjectiveTransform(H).inverse, mode=mode, order=1)
p1 = fast_homography(img, H, mode=mode)
p1 = np.round(p1)
## import matplotlib.pyplot as plt
## f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
## ax0.imshow(img)
## ax1.imshow(p0, cmap=plt.cm.gray)
## ax2.imshow(p1, cmap=plt.cm.gray)
## ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
## plt.show()
# import matplotlib.pyplot as plt
# f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
# ax0.imshow(img)
# ax1.imshow(p0, cmap=plt.cm.gray)
# ax2.imshow(p1, cmap=plt.cm.gray)
# ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
# plt.show()
d = np.mean(np.abs(p0 - p1))
assert d < 0.2
assert d < 0.001
def test_swirl():