I want to speed up my blob detection function by computing the second order derivatives of the input image using the integral image and some box filters.

Initially, I used the following method, which provides correct results.

def smooth_gradient(image, sigma, deg):
    # define the filters according to the arguments
    Gs = myfilter(sigma, "gaussian")
    # smoothen the image
    smooth = cv2.filter2D(image, -1, Gs)
    # calculate the gradient on both directions 
    gradx, grady = np.gradient(smooth)
    if (deg==1):
        return (gradx, grady)
    elif (deg==2):
        gradxx, gradxy = np.gradient(gradx)
        gradxy, gradyy = np.gradient(grady)
        return (gradxx, gradxy, gradyy)
    elif (deg==3):
        gradxx, gradxy = np.gradient(gradx)
        gradxy, gradyy = np.gradient(grady)
        return (gradx, grady, gradxx, gradxy, gradyy)
    print("Usage: deg = 1 for (gradx, grady)")
    print("       deg = 2 for (gradxx, gradxy, gradyy)")
    print("       deg = 3 for (gradx, grady, gradxx, gradxy, gradyy)")

The method using integral images is supposed to be faster, because the traditional method of performing the convolution of the image with some kernel gets slower as the size of the kernel increases, whereas the integral image method is independent of the kernel size. As per my understanding, this happens because the integral image holds all the necessary sums that the convolution would have to constantly recalculate.

I implemented a simple method for computing the integral image of the given image.

def IntegralImage(i):
    return np.cumsum(np.cumsum(i, axis=0), axis=1)

Then, i iterated through the pixels of the integral image. In each iteration, i used the current pixel as the center of a rectangular area and calculated the coordinates of the rectangle's corners.

def BoxDerivative(ii, sigma):
   n = int(2*np.ceil(3*sigma) + 1)
   d1 = int(4*np.floor(n/6) + 1)
   d2 = int(2*np.floor(n/6) + 1)
   (x, y) = ii.shape
   for ix in range(x):
      for iy in range(y):
          # simple operations that take up too many lines ...
          gradxx = (rect area) -2*(another rect area) + (another rect area)
          gradyy = ...
          gradxy = ...

Finally, I calculated the sum of the pixels inside the rectangular area using the values of the integral image at the corners. I did that for some specific rectangles that are supposed to approximate the second order derivatives in each direction (x, y and xy).

I expected this method to be fast, however it is remarkably slower, around 200 times slower than my previous method. I know that the bottleneck is the iteration through the image's pixels. This leads me to believe that my understanding of the integral image and box filter method is lacking.

  • $\begingroup$ Loops in Python are slow. np.gradient does the loops in a compiled language. You need to have a very expensive convolution for the integral image method implemented in pure Python to be faster. $\endgroup$ Apr 2, 2023 at 16:33
  • $\begingroup$ Also, the integral image method is for box filters, not for derivative filters. I don’t think it makes sense to try to use it for those filters. $\endgroup$ Apr 2, 2023 at 16:34
  • $\begingroup$ @CrisLuengo thanks for letting me know that Python loops are slow. Perhaps I was not clear enough. I use box filters that approximate the derivative. Be that as it may, is my understanding of the concept alright? I want the derivative of the image to be an array of the same shape, therefore I believe I have to iterate through the pixels of the integral image and for each pixel compute the box filter that approximates the derivative of each direction, which requires constant time. That’s why I’m looping. $\endgroup$ Apr 2, 2023 at 16:47
  • $\begingroup$ You’re better off using an existing implementation in a library like OpenCV, which is compiled and will run very fast. Write the function in Python to learn how it works if you need to, but don’t expect it to be fast. $\endgroup$ Apr 2, 2023 at 17:23
  • 1
    $\begingroup$ ah, now I see what you're going for, a kind of Haar feature. yes, those can be calculated cheaply from integral images. Haar features can very coarsely approximate various kernel shapes. for small kernels, all this extra trouble isn't worth it. it makes sense for larger kernels, because the lookup in the integral image will cost the same, while direct calculation scales with the area of the kernel. the break-even point, if it matters, would depend on many things and all but requires benchmarks to figure out. $\endgroup$ Jul 8, 2023 at 20:22

1 Answer 1


Summarization of comments into an answer:

  1. Python is interpreted with no JIT acceleration. Loops are slow.
  2. There is a package called Numba which accelerates the work on NumPy arrays (Loops mostly).
  3. Integral Images are fast calculation of the box filter. With Box filter you may approximate Gaussian Filter (See Gaussian Equivalent of Convolving an Image 50 Times with a Box Filter).
  4. There are other accelerated implementations of Box Blur / Box Filter.
  5. Mind numerical issues with Integral Images.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.