diff --git a/bento.info b/bento.info index 7d21bfa8..659247d1 100644 --- a/bento.info +++ b/bento.info @@ -98,6 +98,9 @@ Library: Extension: skimage.feature.corner_cy Sources: skimage/feature/corner_cy.pyx + Extension: skimage.feature._hoghistogram + Sources: + skimage/feature/_hoghistogram.pyx Extension: skimage.feature._texture Sources: skimage/feature/_texture.pyx diff --git a/skimage/feature/_hog.py b/skimage/feature/_hog.py index 11f7099d..afb2b4e9 100644 --- a/skimage/feature/_hog.py +++ b/skimage/feature/_hog.py @@ -1,7 +1,7 @@ +from __future__ import division import numpy as np -from scipy import sqrt, pi, arctan2, cos, sin -from scipy.ndimage import uniform_filter from .._shared.utils import assert_nD +from . import _hoghistogram def hog(image, orientations=9, pixels_per_cell=(8, 8), @@ -63,7 +63,7 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), assert_nD(image, 2) if normalise: - image = sqrt(image) + image = np.sqrt(image) """ The second stage computes first order image gradients. These capture @@ -104,9 +104,6 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), cell are used to vote into the orientation histogram. """ - magnitude = sqrt(gx ** 2 + gy ** 2) - orientation = arctan2(gy, gx) * (180 / pi) % 180 - sy, sx = image.shape cx, cy = pixels_per_cell bx, by = cells_per_block @@ -116,22 +113,9 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), # compute orientations integral images orientation_histogram = np.zeros((n_cellsy, n_cellsx, orientations)) - subsample = np.index_exp[cy // 2:cy * n_cellsy:cy, - cx // 2:cx * n_cellsx:cx] - for i in range(orientations): - # create new integral image for this orientation - # isolate orientations in this range - temp_ori = np.where(orientation < 180.0 / orientations * (i + 1), - orientation, -1) - temp_ori = np.where(orientation >= 180.0 / orientations * i, - temp_ori, -1) - # select magnitudes for those orientations - cond2 = temp_ori > -1 - temp_mag = np.where(cond2, magnitude, 0) - - temp_filt = uniform_filter(temp_mag, size=(cy, cx)) - orientation_histogram[:, :, i] = temp_filt[subsample] + _hoghistogram.hog_histograms(gx, gy, cx, cy, sx, sy, n_cellsx, n_cellsy, + orientations, orientation_histogram) # now for each cell, compute the histogram hog_image = None @@ -140,13 +124,16 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), from .. import draw radius = min(cx, cy) // 2 - 1 + orientations_arr = np.arange(orientations) + dx_arr = radius * np.cos(orientations_arr / orientations * np.pi) + dy_arr = radius * np.sin(orientations_arr / orientations * np.pi) + cr2 = cy + cy + cc2 = cx + cx hog_image = np.zeros((sy, sx), dtype=float) for x in range(n_cellsx): for y in range(n_cellsy): - for o in range(orientations): - centre = tuple([y * cy + cy // 2, x * cx + cx // 2]) - dx = radius * cos(float(o) / orientations * np.pi) - dy = radius * sin(float(o) / orientations * np.pi) + for o, dx, dy in zip(orientations_arr, dx_arr, dy_arr): + centre = tuple([y * cr2 // 2, x * cc2 // 2]) rr, cc = draw.line(int(centre[0] - dx), int(centre[1] + dy), int(centre[0] + dx), @@ -177,7 +164,7 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), for y in range(n_blocksy): block = orientation_histogram[y:y + by, x:x + bx, :] eps = 1e-5 - normalised_blocks[y, x, :] = block / sqrt(block.sum() ** 2 + eps) + normalised_blocks[y, x, :] = block / np.sqrt(block.sum() ** 2 + eps) """ The final step collects the HOG descriptors from all blocks of a dense diff --git a/skimage/feature/_hoghistogram.pyx b/skimage/feature/_hoghistogram.pyx new file mode 100644 index 00000000..3f8899dd --- /dev/null +++ b/skimage/feature/_hoghistogram.pyx @@ -0,0 +1,142 @@ +# cython: cdivision=True +# cython: boundscheck=False +# cython: wraparound=False + +import numpy as np +cimport numpy as cnp + +cdef float cell_hog(double[:, ::1] magnitude, + double[:, ::1] orientation, + float orientation_start, float orientation_end, + int cell_columns, int cell_rows, + int column_index, int row_index, + int size_columns, int size_rows) nogil: + """Calculation of the cell's HOG value + + Parameters + ---------- + magnitude : ndarray + The gradient magnitudes of the pixels. + orientation : ndarray + Lookup table for orientations. + orientation_start : float + Orientation range start. + orientation_end : float + Orientation range end. + cell_columns : int + Pixels per cell (rows). + cell_rows : int + Pixels per cell (columns). + column_index : int + Block column index. + row_index : int + Block row index. + size_columns : int + Number of columns. + size_rows : int + Number of rows. + + Returns + ------- + total : float + The total HOG value. + """ + cdef int cell_column, cell_row, cell_row_index, cell_column_index, \ + range_columns_start, range_columns_stop, range_rows_start, \ + range_rows_stop + + range_rows_stop = cell_rows/2 + range_rows_start = -range_rows_stop + range_columns_stop = cell_columns/2 + range_columns_start = -range_columns_stop + + cdef float total = 0. + for cell_row in range(range_rows_start, range_rows_stop): + cell_row_index = row_index + cell_row + if (cell_row_index < 0 or cell_row_index >= size_rows): + continue + + for cell_column in range(range_columns_start, range_columns_stop): + cell_column_index = column_index + cell_column + if (cell_column_index < 0 or cell_column_index >= size_columns + or orientation[cell_row_index, cell_column_index] + >= orientation_start + or orientation[cell_row_index, cell_column_index] + < orientation_end): + continue + + total += magnitude[cell_row_index, cell_column_index] + + return total + +def hog_histograms(double[:, ::1] gradient_columns, + double[:, ::1] gradient_rows, + int cell_columns, int cell_rows, + int size_columns, int size_rows, + int number_of_cells_columns, int number_of_cells_rows, + int number_of_orientations, + cnp.float64_t[:, :, :] orientation_histogram): + """Extract Histogram of Oriented Gradients (HOG) for a given image. + + Parameters + ---------- + gradient_columns : ndarray + First order image gradients (rows). + gradient_rows : ndarray + First order image gradients (columns). + cell_columns : int + Pixels per cell (rows). + cell_rows : int + Pixels per cell (columns). + size_columns : int + Number of columns. + size_rows : int + Number of rows. + number_of_cells_columns : int + Number of cells (rows). + number_of_cells_rows : int + Number of cells (columns). + number_of_orientations : int + Number of orientation bins. + orientation_histogram : ndarray + The histogram array which is modified in place. + """ + + cdef double[:, ::1] magnitude = np.hypot(gradient_columns, + gradient_rows) + cdef double[:, ::1] orientation = \ + np.arctan2(gradient_rows, gradient_columns) * (180 / np.pi) % 180 + cdef int i, x, y, o, yi, xi, cc, cr, x0, y0 + cdef float orientation_start, orientation_end, \ + number_of_orientations_per_180 + + x0 = cell_columns / 2 + y0 = cell_rows / 2 + cc = cell_rows * number_of_cells_rows + cr = cell_columns * number_of_cells_columns + number_of_orientations_per_180 = 180. / number_of_orientations + + with nogil: + # compute orientations integral images + for i in range(number_of_orientations): + # isolate orientations in this range + orientation_start = number_of_orientations_per_180 * (i + 1) + orientation_end = number_of_orientations_per_180 * i + x = x0 + y = y0 + yi = 0 + xi = 0 + + while y < cc: + xi = 0 + x = x0 + + while x < cr: + orientation_histogram[yi, xi, i] = cell_hog(magnitude, + orientation, orientation_start, orientation_end, + cell_columns, cell_rows, x, y, size_columns, size_rows) + xi += 1 + x += cell_columns + + yi += 1 + y += cell_rows diff --git a/skimage/feature/setup.py b/skimage/feature/setup.py index d7be3dcf..62bd950b 100644 --- a/skimage/feature/setup.py +++ b/skimage/feature/setup.py @@ -18,6 +18,7 @@ def configuration(parent_package='', top_path=None): cython(['brief_cy.pyx'], working_path=base_path) cython(['_texture.pyx'], working_path=base_path) cython(['_hessian_det_appx.pyx'], working_path=base_path) + cython(['_hoghistogram.pyx'], working_path=base_path) config.add_extension('corner_cy', sources=['corner_cy.c'], include_dirs=[get_numpy_include_dirs()]) @@ -31,6 +32,8 @@ def configuration(parent_package='', top_path=None): include_dirs=[get_numpy_include_dirs(), '../_shared']) config.add_extension('_hessian_det_appx', sources=['_hessian_det_appx.c'], include_dirs=[get_numpy_include_dirs()]) + config.add_extension('_hoghistogram', sources=['_hoghistogram.c'], + include_dirs=[get_numpy_include_dirs(), '../_shared']) return config