From 2f897dc35dd9ef1329b4be41cd3dabbebc0689ae Mon Sep 17 00:00:00 2001
From: Brian Holt <bdholt1@gmail.com>
Date: Tue, 13 Sep 2011 15:53:01 +0100
Subject: [PATCH 1/8] Initial implementation of Histograms of Oriented
 Gradients

---
 scikits/image/feature/__init__.py       |   2 +
 scikits/image/feature/hog.py            | 199 ++++++++++++++++++++++++
 scikits/image/feature/tests/test_hog.py |  25 +++
 3 files changed, 226 insertions(+)
 create mode 100644 scikits/image/feature/__init__.py
 create mode 100644 scikits/image/feature/hog.py
 create mode 100644 scikits/image/feature/tests/test_hog.py

diff --git a/scikits/image/feature/__init__.py b/scikits/image/feature/__init__.py
new file mode 100644
index 00000000..47e14c9d
--- /dev/null
+++ b/scikits/image/feature/__init__.py
@@ -0,0 +1,2 @@
+
+from hog import histogram_of_oriented_gradients
\ No newline at end of file
diff --git a/scikits/image/feature/hog.py b/scikits/image/feature/hog.py
new file mode 100644
index 00000000..7cfc0fa1
--- /dev/null
+++ b/scikits/image/feature/hog.py
@@ -0,0 +1,199 @@
+"""Extract Histogram of Oriented Gradients feature from image."""
+
+# Authors: Brian Holt
+#
+# License: BSD
+
+import numpy as np
+from scipy import sqrt, pi, arctan, cos, sin
+
+
+def histogram_of_oriented_gradients(image, n_orientations=9, ppc=(8,8), 
+                                    cpb=(3,3), visualise=False, 
+                                    apply_normalisation=False):   
+    """ 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 
+
+    Parameters
+    ----------
+    image: ndarray, 2D
+        2D image (greyscale)
+
+    n_orientations  : int
+        number of orientation bins
+
+    ppc  : 2 tuple (int,int)
+        pixels per cell, size in pixels of a cell
+
+    cpb  : 2 tuple (int,int)
+        cells per block, number of cells in each block
+
+    visualise : bool
+        return an image of the HOG
+
+    apply_normalisation : bool
+        apply the initial optional global normalisation
+        
+    Returns
+    -------
+    newarr : ndarray 
+        HOG for the image as a 1D (flattened) array.
+
+    References
+    ----------
+    * http://en.wikipedia.org/wiki/Histogram_of_oriented_gradients  
+    
+    * 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)   
+       
+    """
+    The first stage applies an optional global image normalisation 
+    equalisation that is designed to reduce the influence of illumination 
+    effects. In practice we use gamma (power law) compression, either 
+    computing the square root or the log of each colour channel.
+    Image texture strength is typically proportional to the local surface 
+    illumination so this compression helps to reduce the effects of local 
+    shadowing and illumination variations.
+    """
+    
+    if apply_normalisation:  
+        image = sqrt(image)   
+       
+    """
+    The second stage computes first order image gradients. These capture 
+    contour, silhouette and some texture information, while providing 
+    further resistance to illumination variations. The locally dominant 
+    colour channel is used, which provides colour invariance to a large extent.
+    Variant methods may also include second order image derivatives, which
+    act as primitive bar detectors - a useful feature for capturing, 
+    e.g. bar like structures in bicycles and limbs in humans   
+    """    
+       
+    if image.ndim == 3:
+        # replace RGB with locally dominant colour channel
+        pass # TODO   
+    gx = np.zeros(image.shape)
+    gy = np.zeros(image.shape)       
+    gx[:-1, :-1] = image[:-1,:-1]-image[:-1,1:]
+    gy[:-1, :-1] = image[:-1,:-1]-image[1:,:-1]
+
+    #import Image
+    #Image.fromarray(gx).show()
+    #Image.fromarray(gy).show()    
+    
+
+    """
+    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. 
+    """
+    
+    magnitude = sqrt(gx**2+gy**2)
+    orientation = arctan(gy/(gx+1e-15))*(180/pi)+90 
+
+    # compute n_orientations integral images
+    integral_images = []
+    for i in range(0, n_orientations):
+        #create new integral image for this orientation
+        # isolate orientations in this range
+
+        temp_ori = np.where(orientation < 180/n_orientations*(i+1), 
+                            orientation, 0)
+        temp_ori = np.where(orientation >= 180/n_orientations*i, 
+                            temp_ori, 0)
+        # select magnitudes for those orientations
+        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)
+              
+    sx,sy = image.shape
+    cx, cy = ppc        
+    bx, by = cpb        
+                    
+    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, n_orientations))
+        
+    radius = min(cx, cy) // 2 - 1
+    hog_image = None
+    if visualise:
+        import Image, ImageDraw
+        hog_image = Image.new("F", (sy,sx))
+        draw = ImageDraw.Draw(hog_image)
+                 
+    for x in range(0, n_cellsx):
+        for y in range(0, n_cellsy):
+            for o in range(0, n_orientations):
+                # compute the histogram from integral image
+                #print x, y, o
+                A = integral_images[o][x*cx, y*cy]
+                B = integral_images[o][(x+1)*cx-1, y*cy]
+                C = integral_images[o][(x+1)*cx-1, (y+1)*cy-1]
+                D = integral_images[o][x*cx, (y+1)*cy-1]                
+                orientation_histogram[x, y, o] = A + C - D - B
+                            
+                if visualise:                            
+                    centre = tuple([y*cy + cy//2 , x*cx + cx//2])
+                    dx = radius*cos(1.0*o/n_orientations*np.pi)
+                    dy = radius*sin(1.0*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])
+                        
+    """
+    The fourth stage computes normalisation, which takes local groups of 
+    cells and contrast normalises their overall responses before passing 
+    to next stage. Normalisation introduces better invariance to illumination, 
+    shadowing, and edge contrast. It is performed by accumulating a measure 
+    of local histogram "energy" over local groups of cells that we call 
+    "blocks". The result is used to normalise each cell in the block. 
+    Typically each individual cell is shared between several blocks, but 
+    its normalisations are block dependent and thus different. The cell 
+    thus appears several times in the final output vector with different 
+    normalisations. This may seem redundant but it improves the performance. 
+    We refer to the normalised block descriptors as Histogram of Oriented 
+    Gradient (HOG) descriptors.
+    """            
+    
+    n_blocksx = (n_cellsx - bx) + 1
+    n_blocksy = (n_cellsy - by) + 1     
+    normalised_blocks = np.zeros((n_blocksx, n_blocksy, 
+                                  bx, by, n_orientations))
+    
+    for x in range(0, n_blocksx):
+        for y in range(0, 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
+    """
+    
+    return normalised_blocks.ravel(), hog_image    
diff --git a/scikits/image/feature/tests/test_hog.py b/scikits/image/feature/tests/test_hog.py
new file mode 100644
index 00000000..ba3de467
--- /dev/null
+++ b/scikits/image/feature/tests/test_hog.py
@@ -0,0 +1,25 @@
+# Authors: Brian Holt
+#
+# License: BSD
+
+import numpy as np
+import scipy as sp
+from scipy import ndimage
+
+from numpy.testing import assert_raises
+
+from scikits.image.feature import histogram_of_oriented_gradients  
+
+def test_histogram_of_oriented_gradients():
+    img = sp.lena().astype(np.int8) 
+    
+    fd, hog_image = histogram_of_oriented_gradients(img, 
+                                                    n_orientations=9, 
+                                                    ppc=(8,8), 
+                                                    cpb=(1,1), 
+                                                    visualise=False)
+    assert len(fd) == 9 * (512//8) ** 2
+    
+if __name__ == '__main__':
+    import nose
+    nose.runmodule()

From 207f8661f163b300fd01cc29098b6c0ddaee9a84 Mon Sep 17 00:00:00 2001
From: Brian Holt <bdholt1@gmail.com>
Date: Sat, 17 Sep 2011 18:35:03 +0100
Subject: [PATCH 2/8] pep8 + addressing some of Stefan's comments

---
 CONTRIBUTORS.txt                        |   3 +
 scikits/image/feature/__init__.py       |   3 +-
 scikits/image/feature/hog.py            | 231 ++++++++++++------------
 scikits/image/feature/tests/test_hog.py |  24 +--
 4 files changed, 130 insertions(+), 131 deletions(-)

diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index bebd63ed..995bfc26 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -59,3 +59,6 @@
 
 - Kyle Mandli
   CSV to ReST code for MATLAB comparison table.
+
+- Brian Holt
+  Histograms of Oriented Gradients
\ No newline at end of file
diff --git a/scikits/image/feature/__init__.py b/scikits/image/feature/__init__.py
index 47e14c9d..eb6faa62 100644
--- a/scikits/image/feature/__init__.py
+++ b/scikits/image/feature/__init__.py
@@ -1,2 +1 @@
-
-from hog import histogram_of_oriented_gradients
\ No newline at end of file
+from hog import hog
\ No newline at end of file
diff --git a/scikits/image/feature/hog.py b/scikits/image/feature/hog.py
index 7cfc0fa1..8d9dc4f7 100644
--- a/scikits/image/feature/hog.py
+++ b/scikits/image/feature/hog.py
@@ -1,24 +1,22 @@
-"""Extract Histogram of Oriented Gradients feature from image."""
-
-# Authors: Brian Holt
-#
-# License: BSD
+"""
+:author: Brian Holt, 2011
+:license: modified BSD
+"""
 
 import numpy as np
 from scipy import sqrt, pi, arctan, cos, sin
 
 
-def histogram_of_oriented_gradients(image, n_orientations=9, ppc=(8,8), 
-                                    cpb=(3,3), visualise=False, 
-                                    apply_normalisation=False):   
-    """ Histogram of oriented gradients (HOG) for a given image.
+def hog(image, n_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.
 
     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 
+        4) flatten into a feature vector
 
     Parameters
     ----------
@@ -28,86 +26,91 @@ def histogram_of_oriented_gradients(image, n_orientations=9, ppc=(8,8),
     n_orientations  : int
         number of orientation bins
 
-    ppc  : 2 tuple (int,int)
+    pixels_per_cell  : 2 tuple (int,int)
         pixels per cell, size in pixels of a cell
 
-    cpb  : 2 tuple (int,int)
+    cells_per_block  : 2 tuple (int,int)
         cells per block, number of cells in each block
 
-    visualise : bool
+    visualise : bool, optional
         return an image of the HOG
 
-    apply_normalisation : bool
-        apply the initial optional global normalisation
-        
+    normalise : bool, optional
+        apply power law compression to normalise the image before
+        processing
+
     Returns
     -------
-    newarr : ndarray 
+    newarr : ndarray
         HOG for the image as a 1D (flattened) array.
 
+    hog_image : PIL Image (if visualise=True)
+        A visualisation of the HOG image
+
     References
     ----------
-    * http://en.wikipedia.org/wiki/Histogram_of_oriented_gradients  
-    
-    * 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
+    * http://en.wikipedia.org/wiki/Histogram_of_oriented_gradients
+
+    * 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)   
-       
+
+    image = np.atleast_2d(image)
+
     """
-    The first stage applies an optional global image normalisation 
-    equalisation that is designed to reduce the influence of illumination 
-    effects. In practice we use gamma (power law) compression, either 
+    The first stage applies an optional global image normalisation
+    equalisation that is designed to reduce the influence of illumination
+    effects. In practice we use gamma (power law) compression, either
     computing the square root or the log of each colour channel.
-    Image texture strength is typically proportional to the local surface 
-    illumination so this compression helps to reduce the effects of local 
+    Image texture strength is typically proportional to the local surface
+    illumination so this compression helps to reduce the effects of local
     shadowing and illumination variations.
     """
-    
-    if apply_normalisation:  
-        image = sqrt(image)   
-       
-    """
-    The second stage computes first order image gradients. These capture 
-    contour, silhouette and some texture information, while providing 
-    further resistance to illumination variations. The locally dominant 
-    colour channel is used, which provides colour invariance to a large extent.
-    Variant methods may also include second order image derivatives, which
-    act as primitive bar detectors - a useful feature for capturing, 
-    e.g. bar like structures in bicycles and limbs in humans   
-    """    
-       
+
     if image.ndim == 3:
         # replace RGB with locally dominant colour channel
-        pass # TODO   
+        pass  # TODO
+    if normalise:
+        image = sqrt(image)
+
+    """
+    The second stage computes first order image gradients. These capture
+    contour, silhouette and some texture information, while providing
+    further resistance to illumination variations. The locally dominant
+    colour channel is used, which provides colour invariance to a large
+    extent. Variant methods may also include second order image derivatives,
+    which act as primitive bar detectors - a useful feature for capturing,
+    e.g. bar like structures in bicycles and limbs in humans.
+    """
+
     gx = np.zeros(image.shape)
-    gy = np.zeros(image.shape)       
+    gy = np.zeros(image.shape)
     gx[:-1, :-1] = image[:-1,:-1]-image[:-1,1:]
     gy[:-1, :-1] = image[:-1,:-1]-image[1:,:-1]
+    #gx[:-1, :-1] = np.diff(image, n=1, axis=0)
+    #gy[:-1, :-1] = np.diff(image, n=1, axis=1)
 
     #import Image
     #Image.fromarray(gx).show()
-    #Image.fromarray(gy).show()    
-    
+    #Image.fromarray(gy).show()
 
     """
-    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 
+    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. 
+    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 = arctan(gy/(gx+1e-15))*(180/pi)+90 
+
+    magnitude = sqrt(gx ** 2 + gy ** 2)
+    orientation = arctan(gy / (gx + 1e-15)) * (180 / pi) + 90
 
     # compute n_orientations integral images
     integral_images = []
@@ -115,85 +118,89 @@ def histogram_of_oriented_gradients(image, n_orientations=9, ppc=(8,8),
         #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 / n_orientations * (i + 1),
                             orientation, 0)
-        temp_ori = np.where(orientation >= 180/n_orientations*i, 
+        temp_ori = np.where(orientation >= 180 / n_orientations * i,
                             temp_ori, 0)
         # select magnitudes for those orientations
         cond2 = temp_ori > 0
-        temp_mag = np.where(cond2, magnitude, 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)    
+                axis=1, dtype=float)
         integral_images.append(integral)
-              
-    sx,sy = image.shape
-    cx, cy = ppc        
-    bx, by = cpb        
-                    
-    n_cellsx = int(np.floor(sx//cx))  # number of cells in x
-    n_cellsy = int(np.floor(sy//cy))  # number of cells in y
-        
+
+    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, n_orientations))
-        
+
     radius = min(cx, cy) // 2 - 1
     hog_image = None
     if visualise:
-        import Image, ImageDraw
-        hog_image = Image.new("F", (sy,sx))
+        import Image
+        import ImageDraw
+        hog_image = Image.new("F", (sy, sx))
         draw = ImageDraw.Draw(hog_image)
-                 
+
     for x in range(0, n_cellsx):
         for y in range(0, n_cellsy):
             for o in range(0, n_orientations):
                 # compute the histogram from integral image
                 #print x, y, o
-                A = integral_images[o][x*cx, y*cy]
-                B = integral_images[o][(x+1)*cx-1, y*cy]
-                C = integral_images[o][(x+1)*cx-1, (y+1)*cy-1]
-                D = integral_images[o][x*cx, (y+1)*cy-1]                
+                A = integral_images[o][x * cx, y * cy]
+                B = integral_images[o][(x + 1) * cx - 1, y * cy]
+                C = integral_images[o][(x + 1) * cx - 1, (y + 1) * cy - 1]
+                D = integral_images[o][x * cx, (y + 1) * cy - 1]
                 orientation_histogram[x, y, o] = A + C - D - B
-                            
-                if visualise:                            
-                    centre = tuple([y*cy + cy//2 , x*cx + cx//2])
-                    dx = radius*cos(1.0*o/n_orientations*np.pi)
-                    dy = radius*sin(1.0*o/n_orientations*np.pi)
-                    draw.line([(centre[0]-dx, centre[1]-dy), 
-                               (centre[0]+dx, centre[1]+dy)], 
+
+                if visualise:
+                    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])
-                        
+
     """
-    The fourth stage computes normalisation, which takes local groups of 
-    cells and contrast normalises their overall responses before passing 
-    to next stage. Normalisation introduces better invariance to illumination, 
-    shadowing, and edge contrast. It is performed by accumulating a measure 
-    of local histogram "energy" over local groups of cells that we call 
-    "blocks". The result is used to normalise each cell in the block. 
-    Typically each individual cell is shared between several blocks, but 
-    its normalisations are block dependent and thus different. The cell 
-    thus appears several times in the final output vector with different 
-    normalisations. This may seem redundant but it improves the performance. 
-    We refer to the normalised block descriptors as Histogram of Oriented 
+    The fourth stage computes normalisation, which takes local groups of
+    cells and contrast normalises their overall responses before passing
+    to next stage. Normalisation introduces better invariance to illumination,
+    shadowing, and edge contrast. It is performed by accumulating a measure
+    of local histogram "energy" over local groups of cells that we call
+    "blocks". The result is used to normalise each cell in the block.
+    Typically each individual cell is shared between several blocks, but
+    its normalisations are block dependent and thus different. The cell
+    thus appears several times in the final output vector with different
+    normalisations. This may seem redundant but it improves the performance.
+    We refer to the normalised block descriptors as Histogram of Oriented
     Gradient (HOG) descriptors.
-    """            
-    
+    """
+
     n_blocksx = (n_cellsx - bx) + 1
-    n_blocksy = (n_cellsy - by) + 1     
-    normalised_blocks = np.zeros((n_blocksx, n_blocksy, 
+    n_blocksy = (n_cellsy - by) + 1
+    normalised_blocks = np.zeros((n_blocksx, n_blocksy,
                                   bx, by, n_orientations))
-    
+
     for x in range(0, n_blocksx):
         for y in range(0, n_blocksy):
-            block = orientation_histogram[x:x+bx, y:y+by, :]
+            block = orientation_histogram[x:x + bx, y:y + by, :]
             eps = 1e-5
-            normalised_blocks[x, y, :] = block / sqrt(block.sum()**2 + eps) 
-    
+            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 
+     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
     """
-    
-    return normalised_blocks.ravel(), hog_image    
+
+    if visualise:
+        return normalised_blocks.ravel(), hog_image
+    else:
+        return normalised_blocks.ravel()
diff --git a/scikits/image/feature/tests/test_hog.py b/scikits/image/feature/tests/test_hog.py
index ba3de467..b37e0e98 100644
--- a/scikits/image/feature/tests/test_hog.py
+++ b/scikits/image/feature/tests/test_hog.py
@@ -1,25 +1,15 @@
-# Authors: Brian Holt
-#
-# License: BSD
-
 import numpy as np
-import scipy as sp
-from scipy import ndimage
+import scipy
 
-from numpy.testing import assert_raises
-
-from scikits.image.feature import histogram_of_oriented_gradients  
+from scikits.image.feature import hog 
 
 def test_histogram_of_oriented_gradients():
-    img = sp.lena().astype(np.int8) 
+    img = scipy.lena().astype(np.int8) 
     
-    fd, hog_image = histogram_of_oriented_gradients(img, 
-                                                    n_orientations=9, 
-                                                    ppc=(8,8), 
-                                                    cpb=(1,1), 
-                                                    visualise=False)
+    fd = hog(img, n_orientations=9, pixels_per_cell=(8, 8), 
+             cells_per_block=(1, 1))
     assert len(fd) == 9 * (512//8) ** 2
     
 if __name__ == '__main__':
-    import nose
-    nose.runmodule()
+    from numpy.testing import run_module_suite
+    run_module_suite()

From 55c77b54df9db06164fbebde937785fdae68ec01 Mon Sep 17 00:00:00 2001
From: Brian Holt <bdholt1@gmail.com>
Date: Sat, 17 Sep 2011 18:55:41 +0100
Subject: [PATCH 3/8] replaced gradient operation with np.diff

---
 scikits/image/feature/hog.py | 10 ++--------
 1 file changed, 2 insertions(+), 8 deletions(-)

diff --git a/scikits/image/feature/hog.py b/scikits/image/feature/hog.py
index 8d9dc4f7..bdef8d2a 100644
--- a/scikits/image/feature/hog.py
+++ b/scikits/image/feature/hog.py
@@ -86,14 +86,8 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8),
 
     gx = np.zeros(image.shape)
     gy = np.zeros(image.shape)
-    gx[:-1, :-1] = image[:-1,:-1]-image[:-1,1:]
-    gy[:-1, :-1] = image[:-1,:-1]-image[1:,:-1]
-    #gx[:-1, :-1] = np.diff(image, n=1, axis=0)
-    #gy[:-1, :-1] = np.diff(image, n=1, axis=1)
-
-    #import Image
-    #Image.fromarray(gx).show()
-    #Image.fromarray(gy).show()
+    gx[:, :-1] = np.diff(image, n=1, axis=1)
+    gy[:-1, :] = np.diff(image, n=1, axis=0)
 
     """
     The third stage aims to produce an encoding that is sensitive to

From 0d50582ccbfcaf3e9b6ae875f57893de6fbbcb92 Mon Sep 17 00:00:00 2001
From: Brian Holt <bdholt1@gmail.com>
Date: Sat, 17 Sep 2011 18:57:36 +0100
Subject: [PATCH 4/8] replaced arctan with arctan2

---
 scikits/image/feature/hog.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/scikits/image/feature/hog.py b/scikits/image/feature/hog.py
index bdef8d2a..cf4dd852 100644
--- a/scikits/image/feature/hog.py
+++ b/scikits/image/feature/hog.py
@@ -4,7 +4,7 @@
 """
 
 import numpy as np
-from scipy import sqrt, pi, arctan, cos, sin
+from scipy import sqrt, pi, arctan2, cos, sin
 
 
 def hog(image, n_orientations=9, pixels_per_cell=(8, 8),
@@ -104,7 +104,7 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8),
     """
 
     magnitude = sqrt(gx ** 2 + gy ** 2)
-    orientation = arctan(gy / (gx + 1e-15)) * (180 / pi) + 90
+    orientation = arctan2(gy, (gx + 1e-15)) * (180 / pi) + 90
 
     # compute n_orientations integral images
     integral_images = []

From 19b8deb95cd7010af0950b09968a4f84a532652d Mon Sep 17 00:00:00 2001
From: Brian Holt <bdholt1@gmail.com>
Date: Mon, 19 Sep 2011 07:06:12 +0100
Subject: [PATCH 5/8] use sat_sum for performance improvement

---
 scikits/image/feature/hog.py | 13 ++++++-------
 1 file changed, 6 insertions(+), 7 deletions(-)

diff --git a/scikits/image/feature/hog.py b/scikits/image/feature/hog.py
index cf4dd852..8826ddba 100644
--- a/scikits/image/feature/hog.py
+++ b/scikits/image/feature/hog.py
@@ -5,7 +5,7 @@
 
 import numpy as np
 from scipy import sqrt, pi, arctan2, cos, sin
-
+from ..transform import sat_sum
 
 def hog(image, n_orientations=9, pixels_per_cell=(8, 8),
         cells_per_block=(3, 3), visualise=False, normalise=False):
@@ -147,12 +147,11 @@ def hog(image, n_orientations=9, pixels_per_cell=(8, 8),
         for y in range(0, n_cellsy):
             for o in range(0, n_orientations):
                 # compute the histogram from integral image
-                #print x, y, o
-                A = integral_images[o][x * cx, y * cy]
-                B = integral_images[o][(x + 1) * cx - 1, y * cy]
-                C = integral_images[o][(x + 1) * cx - 1, (y + 1) * cy - 1]
-                D = integral_images[o][x * cx, (y + 1) * cy - 1]
-                orientation_histogram[x, y, o] = A + C - D - B
+                orientation_histogram[x, y, o] = sat_sum(integral_images[o],
+                                                         y * cy,
+                                                         x * cx,
+                                                         (y + 1) * cy - 1,
+                                                         (x + 1) * cx - 1)
 
                 if visualise:
                     centre = tuple([y * cy + cy // 2, x * cx + cx // 2])

From 5fc3d856e24eb733c636bf4055da4a1c508594dd Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Mon, 26 Sep 2011 00:40:23 -0700
Subject: [PATCH 6/8] 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 <stefan@sun.ac.za>
Date: Mon, 26 Sep 2011 00:57:55 -0700
Subject: [PATCH 7/8] 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 <andy@marvin>
Date: Tue, 27 Sep 2011 19:31:37 +0200
Subject: [PATCH 8/8] 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