From 5fc3d856e24eb733c636bf4055da4a1c508594dd Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Mon, 26 Sep 2011 00:40:23 -0700 Subject: [PATCH 1/3] Add drawing module with Bresenham line. --- CONTRIBUTORS.txt | 3 ++ TASKS.txt | 4 ++ scikits/image/draw/__init__.py | 1 + scikits/image/draw/_draw.pyx | 66 +++++++++++++++++++++++++++ scikits/image/draw/draw.py | 6 +++ scikits/image/draw/setup.py | 30 ++++++++++++ scikits/image/draw/tests/test_draw.py | 52 +++++++++++++++++++++ 7 files changed, 162 insertions(+) create mode 100644 scikits/image/draw/__init__.py create mode 100644 scikits/image/draw/_draw.pyx create mode 100644 scikits/image/draw/draw.py create mode 100644 scikits/image/draw/setup.py create mode 100644 scikits/image/draw/tests/test_draw.py diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt index d5a71f60..fda2d646 100644 --- a/CONTRIBUTORS.txt +++ b/CONTRIBUTORS.txt @@ -63,3 +63,6 @@ - Brian Holt Histograms of Oriented Gradients + +- David-Warde Farley, Sturla Molden + Bresenheim line drawing, from snippets on numpy-discussion. diff --git a/TASKS.txt b/TASKS.txt index 8c99a40a..1d32b373 100644 --- a/TASKS.txt +++ b/TASKS.txt @@ -122,3 +122,7 @@ Implement Algorithms - Graph cut segmentation - Probabilistic Hough transform +Drawing +``````` +- Wu's algorithm for lines and circles + diff --git a/scikits/image/draw/__init__.py b/scikits/image/draw/__init__.py new file mode 100644 index 00000000..907a1ba7 --- /dev/null +++ b/scikits/image/draw/__init__.py @@ -0,0 +1 @@ +from draw import * diff --git a/scikits/image/draw/_draw.pyx b/scikits/image/draw/_draw.pyx new file mode 100644 index 00000000..a586ae03 --- /dev/null +++ b/scikits/image/draw/_draw.pyx @@ -0,0 +1,66 @@ +import numpy as np +cimport numpy as np +cimport cython + +cdef extern from "math.h": + int abs(int i) + +@cython.boundscheck(False) +@cython.wraparound(False) +def bresenham(int y, int x, int y2, int x2): + """ + Generate line pixel coordinates. + + Parameters + ---------- + y, x : int + Starting position (row, column). + y2, x2 : int + End position (row, column). + + Returns + ------- + rr, cc : (N,) ndarray of int + Indices of pixels that belong to the line. + May be used to directly index into an array, e.g. + ``img[rr, cc] = 1``. + + """ + cdef np.ndarray[np.int32_t, ndim=1, mode="c"] rr, cc + + cdef int steep = 0 + cdef int dx = abs(x2 - x) + cdef int dy = abs(y2 - y) + cdef int sx, sy, d, i + + if (x2 - x) > 0: sx = 1 + else: sx = -1 + if (y2 - y) > 0: sy = 1 + else: sy = -1 + if dy > dx: + steep = 1 + x,y = y,x + dx,dy = dy,dx + sx,sy = sy,sx + d = (2 * dy) - dx + + rr = np.zeros(int(dx) + 1, dtype=np.int32) + cc = np.zeros(int(dx) + 1, dtype=np.int32) + + for i in range(dx): + if steep: + rr[i] = x + cc[i] = y + else: + rr[i] = y + cc[i] = x + while d >= 0: + y = y + sy + d = d - (2 * dx) + x = x + sx + d = d + (2 * dy) + + rr[dx] = y2 + cc[dx] = x2 + + return rr, cc diff --git a/scikits/image/draw/draw.py b/scikits/image/draw/draw.py new file mode 100644 index 00000000..9360623a --- /dev/null +++ b/scikits/image/draw/draw.py @@ -0,0 +1,6 @@ +""" +Methods to draw on arrays. + +""" + +from _draw import bresenham diff --git a/scikits/image/draw/setup.py b/scikits/image/draw/setup.py new file mode 100644 index 00000000..3e0cd1a4 --- /dev/null +++ b/scikits/image/draw/setup.py @@ -0,0 +1,30 @@ +#!/usr/bin/env python + +import os +from scikits.image._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('draw', parent_package, top_path) + config.add_data_dir('tests') + + cython(['_draw.pyx'], working_path=base_path) + + config.add_extension('_draw', sources=['_draw.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 = 'Drawing', + url = 'https://github.com/scikits-image/scikits.image', + license = 'SciPy License (BSD Style)', + **(configuration(top_path='').todict()) + ) diff --git a/scikits/image/draw/tests/test_draw.py b/scikits/image/draw/tests/test_draw.py new file mode 100644 index 00000000..06ebaac0 --- /dev/null +++ b/scikits/image/draw/tests/test_draw.py @@ -0,0 +1,52 @@ +from numpy.testing import assert_array_equal +import numpy as np + +from scikits.image.draw import bresenham + +def test_bresenham_horizontal(): + img = np.zeros((10, 10)) + + rr, cc = bresenham(0, 0, 0, 9) + img[rr, cc] = 1 + + img_ = np.zeros((10, 10)) + img_[0, :] = 1 + + assert_array_equal(img, img_) + +def test_bresenham_vertical(): + img = np.zeros((10, 10)) + + rr, cc = bresenham(0, 0, 9, 0) + img[rr, cc] = 1 + + img_ = np.zeros((10, 10)) + img_[:, 0] = 1 + + assert_array_equal(img, img_) + +def test_reverse(): + img = np.zeros((10, 10)) + + rr, cc = bresenham(0, 9, 0, 0) + img[rr, cc] = 1 + + img_ = np.zeros((10, 10)) + img_[0, :] = 1 + + assert_array_equal(img, img_) + +def test_diag(): + img = np.zeros((5, 5)) + + rr, cc = bresenham(0, 0, 4, 4) + img[rr, cc] = 1 + + img_ = np.eye(5) + + assert_array_equal(img, img_) + + +if __name__ == "__main__": + from numpy.testing import run_module_suite + From 68c61c039b566b8f59e7fb40374f69aaee480a81 Mon Sep 17 00:00:00 2001 From: Stefan van der Walt Date: Mon, 26 Sep 2011 00:57:55 -0700 Subject: [PATCH 2/3] ENH: Use new draw module to construct HOG image. --- scikits/image/feature/hog.py | 126 +++++++++++------------- scikits/image/feature/tests/test_hog.py | 6 +- 2 files changed, 64 insertions(+), 68 deletions(-) diff --git a/scikits/image/feature/hog.py b/scikits/image/feature/hog.py index 8826ddba..8b71f8a3 100644 --- a/scikits/image/feature/hog.py +++ b/scikits/image/feature/hog.py @@ -1,51 +1,42 @@ -""" -:author: Brian Holt, 2011 -:license: modified BSD -""" - import numpy as np from scipy import sqrt, pi, arctan2, cos, sin + +# XXX Replace with integral after merge from ..transform import sat_sum -def hog(image, n_orientations=9, pixels_per_cell=(8, 8), +def hog(image, orientations=9, pixels_per_cell=(8, 8), cells_per_block=(3, 3), visualise=False, normalise=False): - """ Extract Histogram of Oriented Gradients (HOG) for a given image. + """Extract Histogram of Oriented Gradients (HOG) for a given image. Compute a Histogram of Oriented Gradients (HOG) by 1) (optional) global image normalisation 2) computing the gradient image in x and y 3) computing gradient histograms - 3) normalise across blocks - 4) flatten into a feature vector + 3) normalising across blocks + 4) flattening into a feature vector Parameters ---------- - image: ndarray, 2D - 2D image (greyscale) - - n_orientations : int - number of orientation bins - - pixels_per_cell : 2 tuple (int,int) - pixels per cell, size in pixels of a cell - + image : (M, N) ndarray + Input image (greyscale). + orientations : int + Number of orientation bins. + pixels_per_cell : 2 tuple (int, int) + Size (in pixels) of a cell. cells_per_block : 2 tuple (int,int) - cells per block, number of cells in each block - + Number of cells in each block. visualise : bool, optional - return an image of the HOG - + Also return an image of the HOG. normalise : bool, optional - apply power law compression to normalise the image before - processing + Apply power law compression to normalise the image before + processing. Returns ------- newarr : ndarray HOG for the image as a 1D (flattened) array. - - hog_image : PIL Image (if visualise=True) - A visualisation of the HOG image + hog_image : ndarray (if visualise=True) + A visualisation of the HOG image. References ---------- @@ -54,8 +45,8 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8), * Dalal, N and Triggs, B, Histograms of Oriented Gradients for Human Detection, IEEE Computer Society Conference on Computer Vision and Pattern Recognition 2005 San Diego, CA, USA - """ + """ image = np.atleast_2d(image) """ @@ -68,9 +59,9 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8), shadowing and illumination variations. """ - if image.ndim == 3: - # replace RGB with locally dominant colour channel - pass # TODO + if image.ndim > 3: + raise ValueError("Currently only supports grey-level images") + if normalise: image = sqrt(image) @@ -91,30 +82,31 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8), """ The third stage aims to produce an encoding that is sensitive to - local image content while remaining resistant to small changes in pose - or appearance. The adopted method pools gradient orientation information - locally in the same way as the SIFT [Lowe 2004] feature. The image window - is divided into small spatial regions, called "cells". For each cell we - accumulate a local 1-D histogram of gradient or edge orientations over - all the pixels in the cell. This combined cell-level 1-D histogram - forms the basic "orientation histogram" representation. Each orientation - histogram divides the gradient angle range into a fixed number of - predetermined bins. The gradient magnitudes of the pixels in the cell - are used to vote into the orientation histogram. + local image content while remaining resistant to small changes in + pose or appearance. The adopted method pools gradient orientation + information locally in the same way as the SIFT [Lowe 2004] + feature. The image window is divided into small spatial regions, + called "cells". For each cell we accumulate a local 1-D histogram + of gradient or edge orientations over all the pixels in the + cell. This combined cell-level 1-D histogram forms the basic + "orientation histogram" representation. Each orientation histogram + divides the gradient angle range into a fixed number of + predetermined bins. The gradient magnitudes of the pixels in the + cell are used to vote into the orientation histogram. """ magnitude = sqrt(gx ** 2 + gy ** 2) orientation = arctan2(gy, (gx + 1e-15)) * (180 / pi) + 90 - # compute n_orientations integral images + # compute orientations integral images integral_images = [] - for i in range(0, n_orientations): + for i in range(orientations): #create new integral image for this orientation # isolate orientations in this range - temp_ori = np.where(orientation < 180 / n_orientations * (i + 1), + temp_ori = np.where(orientation < 180 / orientations * (i + 1), orientation, 0) - temp_ori = np.where(orientation >= 180 / n_orientations * i, + temp_ori = np.where(orientation >= 180 / orientations * i, temp_ori, 0) # select magnitudes for those orientations cond2 = temp_ori > 0 @@ -122,7 +114,7 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8), #compute integral image integral = np.cumsum(np.cumsum(temp_mag, axis=0, dtype=float), - axis=1, dtype=float) + axis=1, dtype=float) integral_images.append(integral) sx, sy = image.shape @@ -133,19 +125,16 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8), n_cellsy = int(np.floor(sy // cy)) # number of cells in y # now for each cell, compute the histogram - orientation_histogram = np.zeros((n_cellsx, n_cellsy, n_orientations)) + orientation_histogram = np.zeros((n_cellsx, n_cellsy, orientations)) radius = min(cx, cy) // 2 - 1 hog_image = None if visualise: - import Image - import ImageDraw - hog_image = Image.new("F", (sy, sx)) - draw = ImageDraw.Draw(hog_image) + hog_image = np.zeros((sy, sx), dtype=float) - for x in range(0, n_cellsx): - for y in range(0, n_cellsy): - for o in range(0, n_orientations): + for x in range(n_cellsx): + for y in range(n_cellsy): + for o in range(orientations): # compute the histogram from integral image orientation_histogram[x, y, o] = sat_sum(integral_images[o], y * cy, @@ -153,13 +142,18 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8), (y + 1) * cy - 1, (x + 1) * cx - 1) - if visualise: + if visualise: + from scikits.image import draw + + 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) / n_orientations * np.pi) - dy = radius * sin(float(o) / n_orientations * np.pi) - draw.line([(centre[0] - dx, centre[1] - dy), - (centre[0] + dx, centre[1] + dy)], - fill=orientation_histogram[x, y, o]) + dx = radius * cos(float(o) / orientations * np.pi) + dy = radius * sin(float(o) / orientations * np.pi) + rr, cc = draw.bresenham(centre[0] - dx, centre[1] - dy, + centre[0] + dx, centre[1] + dy) + hog_image[rr, cc] += orientation_histogram[x, y, o] """ The fourth stage computes normalisation, which takes local groups of @@ -179,18 +173,18 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8), n_blocksx = (n_cellsx - bx) + 1 n_blocksy = (n_cellsy - by) + 1 normalised_blocks = np.zeros((n_blocksx, n_blocksy, - bx, by, n_orientations)) + bx, by, orientations)) - for x in range(0, n_blocksx): - for y in range(0, n_blocksy): + for x in range(n_blocksx): + for y in range(n_blocksy): block = orientation_histogram[x:x + bx, y:y + by, :] eps = 1e-5 normalised_blocks[x, y, :] = block / sqrt(block.sum() ** 2 + eps) """ - The final step collects the HOG descriptors from all blocks of a dense - overlapping grid of blocks covering the detection window into a combined - feature vector for use in the window classifier + The final step collects the HOG descriptors from all blocks of a dense + overlapping grid of blocks covering the detection window into a combined + feature vector for use in the window classifier. """ if visualise: diff --git a/scikits/image/feature/tests/test_hog.py b/scikits/image/feature/tests/test_hog.py index b37e0e98..6f375cea 100644 --- a/scikits/image/feature/tests/test_hog.py +++ b/scikits/image/feature/tests/test_hog.py @@ -4,10 +4,12 @@ import scipy from scikits.image.feature import hog def test_histogram_of_oriented_gradients(): - img = scipy.lena().astype(np.int8) + # Replace with scikits.image.data.lena() after merge + img = scipy.misc.lena().astype(np.int8) - fd = hog(img, n_orientations=9, pixels_per_cell=(8, 8), + fd = hog(img, orientations=9, pixels_per_cell=(8, 8), cells_per_block=(1, 1)) + assert len(fd) == 9 * (512//8) ** 2 if __name__ == '__main__': From 8f6d290e57ab2fe68c931e96bec3598a1675f88d Mon Sep 17 00:00:00 2001 From: andy Date: Tue, 27 Sep 2011 19:31:37 +0200 Subject: [PATCH 3/3] hog: uniform filter instead of integral image --- scikits/image/feature/hog.py | 33 +++++++++++---------------------- 1 file changed, 11 insertions(+), 22 deletions(-) diff --git a/scikits/image/feature/hog.py b/scikits/image/feature/hog.py index 8b71f8a3..53fafa40 100644 --- a/scikits/image/feature/hog.py +++ b/scikits/image/feature/hog.py @@ -1,5 +1,6 @@ import numpy as np from scipy import sqrt, pi, arctan2, cos, sin +from scipy.ndimage import uniform_filter # XXX Replace with integral after merge from ..transform import sat_sum @@ -98,8 +99,15 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), magnitude = sqrt(gx ** 2 + gy ** 2) orientation = arctan2(gy, (gx + 1e-15)) * (180 / pi) + 90 + sx, sy = image.shape + cx, cy = pixels_per_cell + bx, by = cells_per_block + + n_cellsx = int(np.floor(sx // cx)) # number of cells in x + n_cellsy = int(np.floor(sy // cy)) # number of cells in y + # compute orientations integral images - integral_images = [] + orientation_histogram = np.zeros((n_cellsx, n_cellsy, orientations)) for i in range(orientations): #create new integral image for this orientation # isolate orientations in this range @@ -112,36 +120,17 @@ def hog(image, orientations=9, pixels_per_cell=(8, 8), cond2 = temp_ori > 0 temp_mag = np.where(cond2, magnitude, 0) - #compute integral image - integral = np.cumsum(np.cumsum(temp_mag, axis=0, dtype=float), - axis=1, dtype=float) - integral_images.append(integral) + orientation_histogram[:,:,i] = uniform_filter(temp_mag, size=(cx, cy))[cx/2::cx, cy/2::cy].T - sx, sy = image.shape - cx, cy = pixels_per_cell - bx, by = cells_per_block - - n_cellsx = int(np.floor(sx // cx)) # number of cells in x - n_cellsy = int(np.floor(sy // cy)) # number of cells in y # now for each cell, compute the histogram - orientation_histogram = np.zeros((n_cellsx, n_cellsy, orientations)) + #orientation_histogram = np.zeros((n_cellsx, n_cellsy, orientations)) radius = min(cx, cy) // 2 - 1 hog_image = None if visualise: 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): - # compute the histogram from integral image - orientation_histogram[x, y, o] = sat_sum(integral_images[o], - y * cy, - x * cx, - (y + 1) * cy - 1, - (x + 1) * cx - 1) - if visualise: from scikits.image import draw