mirror of
https://github.com/wassname/scikit-image.git
synced 2026-08-11 11:25:30 +08:00
Merge pull request #384 from andersbll/master
Dense DAISY feature description
This commit is contained in:
Johannes Schönberger
@@ -126,3 +126,6 @@
|
||||
|
||||
- Steven Silvester, Karel Zuiderveld
|
||||
Adaptive Histogram Equalization
|
||||
|
||||
- Anders Boesen Lindbo Larsen
|
||||
Dense DAISY feature description, circle perimeter drawing.
|
||||
|
||||
@@ -0,0 +1,28 @@
|
||||
"""
|
||||
===============================
|
||||
Dense DAISY feature description
|
||||
===============================
|
||||
|
||||
The DAISY local image descriptor is based on gradient orientation histograms
|
||||
similar to the SIFT descriptor. It is formulated in a way that allows for fast
|
||||
dense extraction which is useful for e.g. bag-of-features image
|
||||
representations.
|
||||
|
||||
In this example a limited number of DAISY descriptors are extracted at a large
|
||||
scale for illustrative purposes.
|
||||
"""
|
||||
|
||||
from skimage.feature import daisy
|
||||
from skimage import data
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
|
||||
img = data.camera()
|
||||
descs, descs_img = daisy(img, step=180, radius=58, rings=2, histograms=6,
|
||||
orientations=8, visualize=True)
|
||||
|
||||
plt.axis('off')
|
||||
plt.imshow(descs_img)
|
||||
descs_num = descs.shape[0] * descs.shape[1]
|
||||
plt.title('%i DAISY descriptors extracted:' % descs_num)
|
||||
plt.show()
|
||||
@@ -13,7 +13,7 @@ This example shows how to fill several different shapes:
|
||||
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
from skimage.draw import line, polygon, circle, ellipse
|
||||
from skimage.draw import line, polygon, circle, circle_perimeter, ellipse
|
||||
import numpy as np
|
||||
|
||||
|
||||
@@ -42,5 +42,9 @@ img[rr,cc,:] = (255, 255, 0)
|
||||
rr, cc = ellipse(300, 300, 100, 200, img.shape)
|
||||
img[rr,cc,2] = 255
|
||||
|
||||
# circle
|
||||
rr, cc = circle_perimeter(120, 400, 50)
|
||||
img[rr, cc, :] = (255, 0, 255)
|
||||
|
||||
plt.imshow(img)
|
||||
plt.show()
|
||||
plt.show()
|
||||
|
||||
@@ -1,2 +1,2 @@
|
||||
from ._draw import line, polygon, ellipse, circle
|
||||
from ._draw import line, polygon, ellipse, circle, circle_perimeter, set_color
|
||||
bresenham = line
|
||||
|
||||
@@ -187,3 +187,69 @@ def circle(double cy, double cx, double radius, shape=None):
|
||||
``img[rr, cc] = 1``.
|
||||
"""
|
||||
return ellipse(cy, cx, radius, radius, shape)
|
||||
|
||||
|
||||
def circle_perimeter(int cy, int cx, int radius):
|
||||
"""Generate circle perimeter coordinates.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
cy, cx : int
|
||||
Centre coordinate of circle.
|
||||
radius: int
|
||||
Radius of circle.
|
||||
|
||||
Returns
|
||||
-------
|
||||
rr, cc : (N,) ndarray of int
|
||||
Indices of pixels that belong to the circle perimeter.
|
||||
May be used to directly index into an array, e.g.
|
||||
``img[rr, cc] = 1``.
|
||||
|
||||
"""
|
||||
|
||||
cdef list rr = list()
|
||||
cdef list cc = list()
|
||||
|
||||
cdef int x = 0
|
||||
cdef int y = radius
|
||||
cdef int d = 3 - 2 * radius
|
||||
|
||||
while y >= x:
|
||||
rr.extend([y, -y, y, -y, x, -x, x, -x])
|
||||
cc.extend([x, x, -x, -x, y, y, -y, -y])
|
||||
if d < 0:
|
||||
d += 4 * x + 6
|
||||
else:
|
||||
d += 4 * (x - y) + 10
|
||||
y -= 1
|
||||
x += 1
|
||||
|
||||
return np.array(rr) + cy, np.array(cc) + cx
|
||||
|
||||
|
||||
@cython.boundscheck(False)
|
||||
@cython.wraparound(False)
|
||||
def set_color(img, coords, color):
|
||||
"""Set pixel color in the image at the given coordiantes. Coordinates that
|
||||
exceeed the shape of the image will be ignored.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
img : (M, N, D) ndarray
|
||||
Image
|
||||
coords : ((P,) ndarray, (P,) ndarray)
|
||||
Coordinates of pixels to be colored.
|
||||
color : (D,) ndarray
|
||||
Color to be assigned to coordinates in the image.
|
||||
|
||||
Returns
|
||||
-------
|
||||
img : (M, N, D) ndarray
|
||||
The updated image.
|
||||
"""
|
||||
rr, cc = coords
|
||||
rr_inside = np.logical_and(rr >= 0, rr < img.shape[0])
|
||||
cc_inside = np.logical_and(cc >= 0, cc < img.shape[1])
|
||||
inside = np.logical_and(rr_inside, cc_inside)
|
||||
img[rr[inside], cc[inside]] = color
|
||||
|
||||
@@ -1,7 +1,7 @@
|
||||
from numpy.testing import assert_array_equal
|
||||
import numpy as np
|
||||
|
||||
from skimage.draw import line, polygon, circle, ellipse
|
||||
from skimage.draw import line, polygon, circle, circle_perimeter, ellipse
|
||||
|
||||
|
||||
def test_line_horizontal():
|
||||
@@ -150,6 +150,37 @@ def test_circle():
|
||||
assert_array_equal(img, img_)
|
||||
|
||||
|
||||
def test_circle_perimeter():
|
||||
img = np.zeros((15, 15), 'uint8')
|
||||
rr, cc = circle_perimeter(7, 7, 0)
|
||||
img[rr, cc] = 1
|
||||
assert(np.sum(img) == 1)
|
||||
|
||||
img = np.zeros((17, 15), 'uint8')
|
||||
rr, cc = circle_perimeter(7, 7, 7)
|
||||
img[rr, cc] = 1
|
||||
img_ = np.array(
|
||||
[[0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0],
|
||||
[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
|
||||
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
|
||||
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
|
||||
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
|
||||
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
|
||||
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
|
||||
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
|
||||
[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
|
||||
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
|
||||
[0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0],
|
||||
[0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0],
|
||||
[0, 0, 0, 1, 1, 0, 0, 0, 0, 0, 1, 1, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
|
||||
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]]
|
||||
)
|
||||
assert_array_equal(img, img_)
|
||||
|
||||
|
||||
def test_ellipse():
|
||||
img = np.zeros((15, 15), 'uint8')
|
||||
|
||||
|
||||
@@ -1,3 +1,4 @@
|
||||
from ._daisy import daisy
|
||||
from ._hog import hog
|
||||
from .texture import greycomatrix, greycoprops, local_binary_pattern
|
||||
from .peak import peak_local_max
|
||||
|
||||
@@ -0,0 +1,217 @@
|
||||
import numpy as np
|
||||
from scipy import sqrt, pi, arctan2, cos, sin, exp
|
||||
from scipy.ndimage import gaussian_filter
|
||||
import skimage.color
|
||||
from skimage import img_as_float, draw
|
||||
|
||||
|
||||
def daisy(img, step=4, radius=15, rings=3, histograms=8, orientations=8,
|
||||
normalization='l1', sigmas=None, ring_radii=None, visualize=False):
|
||||
'''Extract DAISY feature descriptors densely for the given image.
|
||||
|
||||
DAISY is a feature descriptor similar to SIFT formulated in a way that
|
||||
allows for fast dense extraction. Typically, this is practical for
|
||||
bag-of-features image representations.
|
||||
|
||||
The implementation follows Tola et al. [1]_ but deviate on the following
|
||||
points:
|
||||
|
||||
* Histogram bin contribution are smoothed with a circular Gaussian
|
||||
window over the tonal range (the angular range).
|
||||
* The sigma values of the spatial Gaussian smoothing in this code do not
|
||||
match the sigma values in the original code by Tola et al. [2]_. In
|
||||
their code, spatial smoothing is applied to both the input image and
|
||||
the center histogram. However, this smoothing is not documented in [1]_
|
||||
and, therefore, it is omitted.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
img : (M, N) array
|
||||
Input image (greyscale).
|
||||
step : int, optional
|
||||
Distance between descriptor sampling points.
|
||||
radius : int, optional
|
||||
Radius (in pixels) of the outermost ring.
|
||||
rings : int, optional
|
||||
Number of rings.
|
||||
histograms : int, optional
|
||||
Number of histograms sampled per ring.
|
||||
orientations : int, optional
|
||||
Number of orientations (bins) per histogram.
|
||||
normalization : [ 'l1' | 'l2' | 'daisy' | 'off' ], optional
|
||||
How to normalize the descriptors
|
||||
|
||||
* 'l1': L1-normalization of each descriptor.
|
||||
* 'l2': L2-normalization of each descriptor.
|
||||
* 'daisy': L2-normalization of individual histograms.
|
||||
* 'off': Disable normalization.
|
||||
|
||||
sigmas : 1D array of float, optional
|
||||
Standard deviation of spatial Gaussian smoothing for the center
|
||||
histogram and for each ring of histograms. The array of sigmas should
|
||||
be sorted from the center and out. I.e. the first sigma value defines
|
||||
the spatial smoothing of the center histogram and the last sigma value
|
||||
defines the spatial smoothing of the outermost ring. Specifying sigmas
|
||||
overrides the following parameter.
|
||||
``rings = len(sigmas)-1``
|
||||
|
||||
ring_radii : 1D array of int, optional
|
||||
Radius (in pixels) for each ring. Specifying ring_radii overrides the
|
||||
following two parameters.
|
||||
| ``rings = len(ring_radii)``
|
||||
| ``radius = ring_radii[-1]``
|
||||
|
||||
If both sigmas and ring_radii are given, they must satisfy the
|
||||
following predicate since no radius is needed for the center
|
||||
histogram.
|
||||
``len(ring_radii) == len(sigmas)+1``
|
||||
visualize : bool, optional
|
||||
Generate a visualization of the DAISY descriptors
|
||||
|
||||
|
||||
Returns
|
||||
-------
|
||||
descs : array
|
||||
Grid of DAISY descriptors for the given image as an array
|
||||
dimensionality (P, Q, R) where
|
||||
| ``P = ceil((M-radius*2)/step)``
|
||||
| ``Q = ceil((N-radius*2)/step)``
|
||||
| ``R = (rings*histograms + 1)*orientations``
|
||||
descs_img : (M, N, 3) array (only if visualize==True)
|
||||
Visualization of the DAISY descriptors.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] Tola et al. "Daisy: An efficient dense descriptor applied to wide-
|
||||
baseline stereo." Pattern Analysis and Machine Intelligence, IEEE
|
||||
Transactions on 32.5 (2010): 815-830.
|
||||
.. [2] http://cvlab.epfl.ch/alumni/tola/daisy.html
|
||||
'''
|
||||
|
||||
# Validate image format.
|
||||
if img.ndim > 2:
|
||||
raise ValueError('Only grey-level images are supported.')
|
||||
if img.dtype.kind != 'f':
|
||||
img = img_as_float(img)
|
||||
|
||||
# Validate parameters.
|
||||
if sigmas is not None and ring_radii is not None \
|
||||
and len(sigmas) - 1 != len(ring_radii):
|
||||
raise ValueError('len(sigmas)-1 != len(ring_radii)')
|
||||
if ring_radii is not None:
|
||||
rings = len(ring_radii)
|
||||
radius = ring_radii[-1]
|
||||
if sigmas is not None:
|
||||
rings = len(sigmas) - 1
|
||||
if sigmas is None:
|
||||
sigmas = [radius * (i + 1) / float(2 * rings) for i in range(rings)]
|
||||
if ring_radii is None:
|
||||
ring_radii = [radius * (i + 1) / float(rings) for i in range(rings)]
|
||||
if normalization not in ['l1', 'l2', 'daisy', 'off']:
|
||||
raise ValueError('Invalid normalization method.')
|
||||
|
||||
# Compute image derivatives.
|
||||
dx = np.zeros(img.shape)
|
||||
dy = np.zeros(img.shape)
|
||||
dx[:, :-1] = np.diff(img, n=1, axis=1)
|
||||
dy[:-1, :] = np.diff(img, n=1, axis=0)
|
||||
|
||||
# Compute gradient orientation and magnitude and their contribution
|
||||
# to the histograms.
|
||||
grad_mag = sqrt(dx ** 2 + dy ** 2)
|
||||
grad_ori = arctan2(dy, dx)
|
||||
orientation_kappa = orientations / pi
|
||||
orientation_angles = [2 * o * pi / orientations - pi
|
||||
for o in range(orientations)]
|
||||
hist = np.empty((orientations,) + img.shape, dtype=float)
|
||||
for i, o in enumerate(orientation_angles):
|
||||
# Weigh bin contribution by the circular normal distribution
|
||||
hist[i, :, :] = exp(orientation_kappa * cos(grad_ori - o))
|
||||
# Weigh bin contribution by the gradient magnitude
|
||||
hist[i, :, :] = np.multiply(hist[i, :, :], grad_mag)
|
||||
|
||||
# Smooth orientation histograms for the center and all rings.
|
||||
sigmas = [sigmas[0]] + sigmas
|
||||
hist_smooth = np.empty((rings + 1,) + hist.shape, dtype=float)
|
||||
for i in range(rings + 1):
|
||||
for j in range(orientations):
|
||||
hist_smooth[i, j, :, :] = gaussian_filter(hist[j, :, :],
|
||||
sigma=sigmas[i])
|
||||
|
||||
# Assemble descriptor grid.
|
||||
theta = [2 * pi * j / histograms for j in range(histograms)]
|
||||
desc_dims = (rings * histograms + 1) * orientations
|
||||
descs = np.empty((desc_dims, img.shape[0] - 2 * radius,
|
||||
img.shape[1] - 2 * radius))
|
||||
descs[:orientations, :, :] = hist_smooth[0, :, radius:-radius,
|
||||
radius:-radius]
|
||||
idx = orientations
|
||||
for i in range(rings):
|
||||
for j in range(histograms):
|
||||
y_min = radius + int(round(ring_radii[i] * sin(theta[j])))
|
||||
y_max = descs.shape[1] + y_min
|
||||
x_min = radius + int(round(ring_radii[i] * cos(theta[j])))
|
||||
x_max = descs.shape[2] + x_min
|
||||
descs[idx:idx + orientations, :, :] = hist_smooth[i + 1, :,
|
||||
y_min:y_max,
|
||||
x_min:x_max]
|
||||
idx += orientations
|
||||
descs = descs[:, ::step, ::step]
|
||||
descs = descs.swapaxes(0, 1).swapaxes(1, 2)
|
||||
|
||||
# Normalize descriptors.
|
||||
if normalization != 'off':
|
||||
descs += 1e-10
|
||||
if normalization == 'l1':
|
||||
descs /= np.sum(descs, axis=2)[:, :, np.newaxis]
|
||||
elif normalization == 'l2':
|
||||
descs /= sqrt(np.sum(descs ** 2, axis=2))[:, :, np.newaxis]
|
||||
elif normalization == 'daisy':
|
||||
for i in range(0, desc_dims, orientations):
|
||||
norms = sqrt(np.sum(descs[:, :, i:i + orientations] ** 2,
|
||||
axis=2))
|
||||
descs[:, :, i:i + orientations] /= norms[:, :, np.newaxis]
|
||||
|
||||
if visualize:
|
||||
descs_img = skimage.color.gray2rgb(img)
|
||||
for i in range(descs.shape[0]):
|
||||
for j in range(descs.shape[1]):
|
||||
# Draw center histogram sigma
|
||||
color = (1, 0, 0)
|
||||
desc_y = i * step + radius
|
||||
desc_x = j * step + radius
|
||||
coords = draw.circle_perimeter(desc_y, desc_x, sigmas[0])
|
||||
draw.set_color(descs_img, coords, color)
|
||||
max_bin = np.max(descs[i, j, :])
|
||||
for o_num, o in enumerate(orientation_angles):
|
||||
# Draw center histogram bins
|
||||
bin_size = descs[i, j, o_num] / max_bin
|
||||
dy = sigmas[0] * bin_size * sin(o)
|
||||
dx = sigmas[0] * bin_size * cos(o)
|
||||
coords = draw.line(desc_y, desc_x, desc_y + dy,
|
||||
desc_x + dx)
|
||||
draw.set_color(descs_img, coords, color)
|
||||
for r_num, r in enumerate(ring_radii):
|
||||
color_offset = float(1 + r_num) / rings
|
||||
color = (1 - color_offset, 1, color_offset)
|
||||
for t_num, t in enumerate(theta):
|
||||
# Draw ring histogram sigmas
|
||||
hist_y = desc_y + int(round(r * sin(t)))
|
||||
hist_x = desc_x + int(round(r * cos(t)))
|
||||
coords = draw.circle_perimeter(hist_y, hist_x,
|
||||
sigmas[r_num + 1])
|
||||
draw.set_color(descs_img, coords, color)
|
||||
for o_num, o in enumerate(orientation_angles):
|
||||
# Draw histogram bins
|
||||
bin_size = descs[i, j, orientations + r_num *
|
||||
histograms * orientations +
|
||||
t_num * orientations + o_num]
|
||||
bin_size /= max_bin
|
||||
dy = sigmas[r_num + 1] * bin_size * sin(o)
|
||||
dx = sigmas[r_num + 1] * bin_size * cos(o)
|
||||
coords = draw.line(hist_y, hist_x, hist_y + dy,
|
||||
hist_x + dx)
|
||||
draw.set_color(descs_img, coords, color)
|
||||
return descs, descs_img
|
||||
else:
|
||||
return descs
|
||||
@@ -0,0 +1,95 @@
|
||||
import numpy as np
|
||||
from numpy.testing import assert_raises, assert_almost_equal
|
||||
from numpy import sqrt, ceil
|
||||
|
||||
from skimage import data
|
||||
from skimage import img_as_float
|
||||
from skimage.feature import daisy
|
||||
|
||||
|
||||
def test_daisy_color_image_unsupported_error():
|
||||
img = np.zeros((20, 20, 3))
|
||||
assert_raises(ValueError, daisy, img)
|
||||
|
||||
|
||||
def test_daisy_desc_dims():
|
||||
img = img_as_float(data.lena()[:128, :128].mean(axis=2))
|
||||
rings = 2
|
||||
histograms = 4
|
||||
orientations = 3
|
||||
descs = daisy(img, rings=rings, histograms=histograms,
|
||||
orientations=orientations)
|
||||
assert(descs.shape[2] == (rings * histograms + 1) * orientations)
|
||||
|
||||
rings = 4
|
||||
histograms = 5
|
||||
orientations = 13
|
||||
descs = daisy(img, rings=rings, histograms=histograms,
|
||||
orientations=orientations)
|
||||
assert(descs.shape[2] == (rings * histograms + 1) * orientations)
|
||||
|
||||
|
||||
def test_descs_shape():
|
||||
img = img_as_float(data.lena()[:256, :256].mean(axis=2))
|
||||
radius = 20
|
||||
step = 8
|
||||
descs = daisy(img, radius=radius, step=step)
|
||||
assert(descs.shape[0] == ceil((img.shape[0] - radius * 2) / float(step)))
|
||||
assert(descs.shape[1] == ceil((img.shape[1] - radius * 2) / float(step)))
|
||||
|
||||
img = img[:-1, :-2]
|
||||
radius = 5
|
||||
step = 3
|
||||
descs = daisy(img, radius=radius, step=step)
|
||||
assert(descs.shape[0] == ceil((img.shape[0] - radius * 2) / float(step)))
|
||||
assert(descs.shape[1] == ceil((img.shape[1] - radius * 2) / float(step)))
|
||||
|
||||
|
||||
def test_daisy_incompatible_sigmas_and_radii():
|
||||
img = img_as_float(data.lena()[:128, :128].mean(axis=2))
|
||||
sigmas = [1, 2]
|
||||
radii = [1, 2]
|
||||
assert_raises(ValueError, daisy, img, sigmas=sigmas, ring_radii=radii)
|
||||
|
||||
|
||||
def test_daisy_normalization():
|
||||
img = img_as_float(data.lena()[:64, :64].mean(axis=2))
|
||||
|
||||
descs = daisy(img, normalization='l1')
|
||||
for i in range(descs.shape[0]):
|
||||
for j in range(descs.shape[1]):
|
||||
assert_almost_equal(np.sum(descs[i, j, :]), 1)
|
||||
descs_ = daisy(img)
|
||||
assert_almost_equal(descs, descs_)
|
||||
|
||||
descs = daisy(img, normalization='l2')
|
||||
for i in range(descs.shape[0]):
|
||||
for j in range(descs.shape[1]):
|
||||
assert_almost_equal(sqrt(np.sum(descs[i, j, :] ** 2)), 1)
|
||||
|
||||
orientations = 8
|
||||
descs = daisy(img, orientations=orientations, normalization='daisy')
|
||||
desc_dims = descs.shape[2]
|
||||
for i in range(descs.shape[0]):
|
||||
for j in range(descs.shape[1]):
|
||||
for k in range(0, desc_dims, orientations):
|
||||
assert_almost_equal(sqrt(np.sum(
|
||||
descs[i, j, k:k + orientations] ** 2)), 1)
|
||||
|
||||
img = np.zeros((50, 50))
|
||||
descs = daisy(img, normalization='off')
|
||||
for i in range(descs.shape[0]):
|
||||
for j in range(descs.shape[1]):
|
||||
assert_almost_equal(np.sum(descs[i, j, :]), 0)
|
||||
|
||||
assert_raises(ValueError, daisy, img, normalization='does_not_exist')
|
||||
|
||||
|
||||
def test_daisy_visualization():
|
||||
img = img_as_float(data.lena()[:128, :128].mean(axis=2))
|
||||
descs, descs_img = daisy(img, visualize=True)
|
||||
assert(descs_img.shape == (128, 128, 3))
|
||||
|
||||
if __name__ == '__main__':
|
||||
from numpy import testing
|
||||
testing.run_module_suite()
|
||||
Reference in New Issue
Block a user