"Histogram of Oriented Gradients" (HOG) feature detector for computer vision
As you indicated in the question, you need to vectorize the for loop:

it = product(xrange(n_divs), xrange(n_divs), xrange(cellx), xrange(celly))
for m, n, i, j in it:
    # grad value
    grad = magnit[m * cellx + i, n * celly + j][0]
    # normalized grad value
    norm_grad = grad / img_area
    # Orientation Angle
    angle = angles[m*cellx + i, n*celly+j][0]
    # (-pi,pi) to (0, 2*pi)
    if angle &lt; 0:
        angle += 2 * pi
    nth_bin = floor(float(angle/BIN_RANGE))
    HOG[((m * n_divs + n) * n_bins + int(nth_bin))] += norm_grad


If you look at what this is doing, you're effectively labelling every pixel in the magnit array with a number below n_HOG, and then summing the normalized values for the pixels with each label.

Operations on labelled regions of images are jobs for the scipy.ndimage.measurements module. Here we can use scipy.ndimage.measurements.sum:

bins = (angles[...,0] % (2 * pi) / BIN_RANGE).astype(int)
x, y = np.mgrid[:width, :height]
x = x * n_divs // width
y = y * n_divs // height
labels = (x * n_divs + y) * n_bins + bins
index = np.arange(n_HOG)
HOG = scipy.ndimage.measurements.sum(magnit[...,0], labels, index)
return HOG / img_area


Notes:


I've used % (2 * pi) to get the angles in the range [0, 2π). An alternative that's more like your code would be angles[angles &lt; 0] += 2 * pi but using the modulus is shorter and, I think, clearer.
I postponed the division by img_area until after the summation, because it looks to me as though in the common case n_HOG is much less than img_area and so it's cheaper to do the division later when there are fewer items. (This means that the results differ very slightly from your code, so bear that in mind when checking.)
I measure the vectorized version as being about 60 times faster than your for loop, but it's still not going to be fast enough to run at 30 fps!
I've written angles[...,0] and magnit[...,0] here in order to drop the third axis. But I think it would make more sense if you dropped this axis earlier, before computing angles and magnit, by writing Ix = Ix[...,0] or just Ix = Ix.reshape((height, width)) if you know that the last axis has length 1.


Update

Based on comments, it looks as if you are using Python 2.7, where the division operator / takes the floor of the result if both arguments are integers. So I've changed the code above to use:

x = x * n_divs // width
y = y * n_divs // height


which is portable between Python 2 and Python 3, and simpler than my first attempt:

x = (x / width * n_divs).astype(int)
y = (y / height * n_divs).astype(int)