Sobel filter kernel of large size

Question

I am using a sobel filter of size 3x3 to calculate the image derivative. Looking at some articles on the internet, it seems that kernels for sobel filter for size 5x5 and 7x7 are also common, but I am not able to find their kernel values.

Could someone please let me know the kernel values for sobel filter of size 5x5 and 7x7? Also, if someone could share a method to generate the kernel values, that will be much useful.

Thanks in advance.

Solution 1

UPDATE 23-Apr-2018: it seems that the kernels defined in the link below are not true Sobel kernels (for 5x5 and above) - they may do a reasonable job of edge detection, but they should not be called Sobel kernels. See Daniel’s answer for a more accurate and comprehensive summary. (I will leave this answer here as (a) it is linked to from various places and (b) accepted answers can not easily be deleted.)

Google seems to turn up plenty of results, e.g. http://rsbweb.nih.gov/nih-image/download/user-macros/slowsobel.macro suggests the following kernels for 3x3, 5x5, 7x7 and 9x9:

3x3:

1   0   -1
2   0   -2
1   0   -1

5x5:

2   1   0   -1  -2
3   2   0   -2  -3
4   3   0   -3  -4
3   2   0   -2  -3
2   1   0   -1  -2

7x7:

3   2   1   0   -1  -2  -3
4   3   2   0   -2  -3  -4
5   4   3   0   -3  -4  -5
6   5   4   0   -4  -5  -6
5   4   3   0   -3  -4  -5
4   3   2   0   -2  -3  -4
3   2   1   0   -1  -2  -3

9x9:

4   3   2   1   0   -1  -2  -3  -4
5   4   3   2   0   -2  -3  -4  -5
6   5   4   3   0   -3  -4  -5  -6
7   6   5   4   0   -4  -5  -6  -7
8   7   6   5   0   -5  -6  -7  -8
7   6   5   4   0   -4  -5  -6  -7
6   5   4   3   0   -3  -4  -5  -6
5   4   3   2   0   -2  -3  -4  -5
4   3   2   1   0   -1  -2  -3  -4

Solution 2

Complete solution for arbitrary Sobel kernel sizes and angles

tl;dr: skip down to section 'Examples'

To add another solution, expanding on this document (it's not particularly high quality, but it shows some usable graphics and matrices starting at the bottom of page 2).

Goal

What we're trying to do is estimate the local gradient of the image at position (x,y). The gradient is a vector made up of the components in x and y direction, gx and gy.

Now, imagine we want to approximate the gradient based on our pixel (x,y) and its neighbours as a kernel operation (3x3, 5x5, or whatever size).

Solution idea

We can approximate the gradient by summing over the projections of all neighbor-center pairs onto the gradient direction. (Sobel's kernel is just a particular method of weighting the different contributions, and so is Prewitt, basically).

Explicit intermediate steps for 3x3

This is the local image, central pixel (x,y) marked as 'o' (center)

a b c
d o f
g h i

Let's say we want the gradient in positive x direction. The unit vector in positive x-direction is (1,0) [I'll later use the convention that the positive y direction is DOWN, i.e. (0,1), and that (0,0) is top left of image).]

The vector from o to f ('of' for short) is (1,0). The gradient in direction 'of' is (f - o) / 1 (value of image at pixel here denoted f minus value at center o, divided by distance between those pixels). If we project the unit vector of that particular neighbor gradient onto our desired gradient direction (1,0) via a dot product we get 1. Here is a little table with the contributions of all neighbors, starting with the easier cases. Note that for diagonals, their distance is sqrt2, and the unit vectors in the diagonal directions are 1/sqrt2 * (+/-1, +/-1)

f:   (f-o)/1     * 1
d:   (d-o)/1     * -1       because (-1, 0) dot (1, 0) = -1
b:   (b-o)/1     * 0        because (0, -1) dot (1, 0) = 0
h:   (h-o)/1     * 0        (as per b)
a:   (a-o)/sqrt2 * -1/sqrt2 distance is sqrt2, and 1/sqrt2*(-1,-1) dot (1,0) = -1/sqrt2
c:   (c-o)/sqrt2 * +1/sqrt2   ...
g:   (g-o)/sqrt2 * -1/sqrt2   ...
i:   (i-o)/sqrt2 * +1/sqrt2   ...

edit for clarification: There are two factors of 1/sqrt(2) for the following reason:

1. We are interested in the contribution to the gradient in a specific direction (here x), so we need to project the directional gradient from the center pixel to the neighbor pixel onto the direction we are interested in. This is accomplished by taking the scalar product of the unit vectors in the respective directions, which introduces the first factor 1/L (here 1/sqrt(2) for the diagonals).

2. The gradient measures the infinitesimal change at a point, which we approximate by finite differences. In terms of a linear equation, m = (y2-y1)/(x2-x1). For this reason, the value difference from the center pixel to the neighbor pixel (y2-y1) has to be distributed over their distance (corresponds to x2-x1) in order to get the ascent units per distance unit. This yields a second factor of 1/L (here 1/sqrt(2) for the diagonals)

Ok, now we know the contributions. Let's simplify this expression by combining opposing pairs of pixel contributions. I'll start with d and f:

{(f-o)/1 * 1} + {(d-o)/1 * -1}
= f - o - (d - o)
= f - d

Now the first diagonal:

{(c-o)/sqrt2 * 1/sqrt2} + {(g-o)/sqrt2 * -1/sqrt2}
= (c - o)/2 - (g - o)/2
= (c - g)/2

The second diagonal contributes (i - a)/2. The perpendicular direction contributes zero. Note that all contributions from the central pixel 'o' vanish.

We have now calculated the contributions of all closest neighbours to the gradient in positive x-direction at pixel (x,y), so our total approximation of the gradient in x-direction is simply their sum:

gx(x,y) = f - d + (c - g)/2 + (i - a)/2

We can obtain the same result by using a convolution kernel where the coefficients are written in the place of the corresponding neighbor pixel:

-1/2  0  1/2
 -1   0   1
-1/2  0  1/2

If you don't want to deal with fractions, you multiply this by 2 and get the well-known Sobel 3x3 kernel.

      -1 0 1
G_x = -2 0 2
      -1 0 1

The multiplication by two only serves to get convenient integers. The scaling of your output image is basically arbitrary, most of the time you normalize it to your image range, anyway (to get clearly visible results).

By the same reasoning as above, you get the kernel for the vertical gradient gy by projecting the neighbor contributions onto the unit vector in positive y direction (0,1)

      -1 -2 -1
G_y =  0  0  0
       1  2  1

Formula for kernels of arbitrary size

If you want 5x5 or larger kernels, you only need to pay attention to the distances, e.g.

A B 2 B A
B C 1 C B
2 1 - 1 2
B C 1 C B
A B 2 B A

where

A = 2 * sqrt2
B = sqrt5
C = sqrt2.

If the length of the vector connecting any two pixels is L, the unit vector in that direction has a prefactor of 1/L. For this reason, the contributions of any pixel 'k' to (say) the x-gradient (1,0) can be simplified to "(value difference over squared distance) times (DotProduct of unnormalized direction vector 'ok' with gradient vector, e.g. (1,0) )"

gx_k = (k - o)/(pixel distance^2) ['ok' dot (1,0)].

Because the dot product of the connecting vector with the x unit vector selects the corresponding vector entry, the corresponding G_x kernel entry at position k is just

i / (i*i + j*j)

where i and j are the number of steps from the center pixel to the pixel k in x and y direction. In the above 3x3 calculation, the pixel 'a' would have i = -1 (1 to the left), j = -1 (1 to the top) and hence the 'a' kernel entry is -1 / (1 + 1) = -1/2.

The entries for the G_y kernel are

j/(i*i + j*j). 

If I want integer values for my kernel, I follow these steps:

• check the available range of the output image
• compute highest possible result from applying floating point kernel (i.e. assume max input value under all positive kernel entries, so output value is (sum over all positive kernel values) * (max possible input image value). If you have signed input, you need to consider the negative values as well. Worst case is then the sum of all positive values + sum of all abs values of negative entries (if max input under positives, -max input under negatives). edit: the sum of all abs values has also been aptly called the weight of the kernel
• calculate maximum allowed up-scaling for kernel (without overflowing range of output image)
• for all integer multiples (from 2 to above maximum) of floating point kernel: check which has the lowest sum of absolute round-off errors and use this kernel

So in summary:

Gx_ij = i / (i*i + j*j)
Gy_ij = j / (i*i + j*j)

where i,j is position in the kernel counted from the center. Scale kernel entries as needed to obtain integer numbers (or at least close approximations).

These formulae hold for all kernel sizes.

Examples

          -2/8 -1/5  0  1/5  2/8           -5  -4  0   4   5
          -2/5 -1/2  0  1/2  2/5           -8 -10  0  10   8
G_x (5x5) -2/4 -1/1  0  1/1  2/4  (*20) = -10 -20  0  20  10
          -2/5 -1/2  0  1/2  2/5           -8 -10  0  10   8
          -2/8 -1/5  0  1/5  2/8           -5  -4  0   4   5

Note that the central 3x3 pixels of the 5x5 kernel in float notation are just the 3x3 kernel, i.e. larger kernels represent a continued approximation with additional but lower-weighted data. This continues on to larger kernel sizes:

           -3/18 -2/13 -1/10 0  1/10 2/13 3/18
           -3/13 -2/8  -1/5  0  1/5  2/8  3/13
           -3/10 -2/5  -1/2  0  1/2  2/5  3/10
G_x (7x7)  -3/9  -2/4  -1/1  0  1/1  2/4  3/9 
           -3/10 -2/5  -1/2  0  1/2  2/5  3/10
           -3/13 -2/8  -1/5  0  1/5  2/8  3/13
           -3/18 -2/13 -1/10 0  1/10 2/13 3/18

Exact integer representations become impractical at this point.

As far as I can tell (don't have access to the original paper), the "Sobel" part to this is properly weighting the contributions. The Prewitt solution can be obtained by leaving out the distance weighting and just entering i and j in the kernel as appropriate.

Bonus: Sobel Kernels for arbitrary directions

So we can approximate the x and y components of the image gradient (which is actually a vector, as stated in the very beginning). The gradient in any arbitrary direction alpha (measured mathematically positive, in this case clockwise since positive y is downward) can be obtained by projecting the gradient vector onto the alpha-gradient unit vector.

The alpha-unit vector is (cos alpha, sin alpha). For alpha = 0° you can obtain the result for gx, for alpha = 90° you get gy.

g_alpha = (alpha-unit vector) dot (gx, gy)
        = (cos a, sin a) dot (gx, gy)
        = cos a * gx + sin a * gy

If you bother to write down gx and gy as sums of neighbor contributions, you realize that you can group the resulting long expression by terms that apply to the same neighbor pixel, and then rewrite this as a single convolution kernel with entries

G_alpha_ij = (i * cos a + j * sin a)/(i*i + j*j)

If you want the closest integer approximation, follow the steps outlined above.

Solution 3

Other sources seem to give different definitions of the larger kernels. The Intel IPP library, for example, gives the 5x5 kernel as

1  2 0  -2 -1
4  8 0  -8 -4
6 12 0 -12 -6
4  8 0  -8 -4
1  2 0  -2 -1

Intuitively, this makes more sense to me because you're paying more attention to the elements closer to the centre. It also has a natural definition in terms of the 3x3 kernel which is easy to extend to generate larger kernels. That said, in my brief search I've found 3 different definitions of the 5x5 kernel - so I suspect that (as Paul says) the larger kernels are ad hoc, and so this is by no means the definitive answer.

The 3x3 kernel is the outer product of a smoothing kernel and a gradient kernel, in Matlab this is something like

sob3x3 = [ 1 2 1 ]' * [1 0 -1]

the larger kernels can be defined by convolving the 3x3 kernel with another smoothing kernel

sob5x5 = conv2( [ 1 2 1 ]' * [1 2 1], sob3x3 )

you can repeat the process to get progressively larger kernels

sob7x7 = conv2( [ 1 2 1 ]' * [1 2 1], sob5x5 )
sob9x9 = conv2( [ 1 2 1 ]' * [1 2 1], sob7x7 )
...

there are a lot of other ways of writing it, but I think this explains exactly what is happening best. Basically, you start off with a smoothing kernel in one direction and a finite differences estimate of the derivative in the other and then just apply smoothing until you get the kernel size you want.

Because it's just a series of convolutions, all the nice properties hold, (commutativity, associativity and so forth) which might be useful for your implementation. For example, you can trivially separate the 5x5 kernel into its smoothing and derivative components:

sob5x5 = conv([1 2 1],[1 2 1])' * conv([1 2 1],[-1 0 1])

Note that in order to be a "proper" derivative estimator, the 3x3 Sobel should be scaled by a factor of 1/8:

sob3x3 = 1/8 * [ 1 2 1 ]' * [1 0 -1]

and each larger kernel needs to be scaled by an additional factor of 1/16 (because the smoothing kernels are not normalised):

sob5x5 = 1/16 * conv2( [ 1 2 1 ]' * [1 2 1], sob3x3 )
sob7x7 = 1/16 * conv2( [ 1 2 1 ]' * [1 2 1], sob5x5 )
...

Solution 4

Here is a simple solution made with python 3 using numpy and the @Daniel answer.

def custom_sobel(shape, axis):
    """
    shape must be odd: eg. (5,5)
    axis is the direction, with 0 to positive x and 1 to positive y
    """
    k = np.zeros(shape)
    p = [(j,i) for j in range(shape[0]) 
           for i in range(shape[1]) 
           if not (i == (shape[1] -1)/2. and j == (shape[0] -1)/2.)]

    for j, i in p:
        j_ = int(j - (shape[0] -1)/2.)
        i_ = int(i - (shape[1] -1)/2.)
        k[j,i] = (i_ if axis==0 else j_)/float(i_*i_ + j_*j_)
    return k

It returns the kernel (5,5) like this:

Sobel x:
   [[-0.25 -0.2   0.    0.2   0.25]
    [-0.4  -0.5   0.    0.5   0.4 ]
    [-0.5  -1.    0.    1.    0.5 ]
    [-0.4  -0.5   0.    0.5   0.4 ]
    [-0.25 -0.2   0.    0.2   0.25]]


Sobel y:
   [[-0.25 -0.4  -0.5  -0.4  -0.25]
    [-0.2  -0.5  -1.   -0.5  -0.2 ]
    [ 0.    0.    0.    0.    0.  ]
    [ 0.2   0.5   1.    0.5   0.2 ]
    [ 0.25  0.4   0.5   0.4   0.25]]

If anyone know a better way to do that in python, please let me know. I'm a newbie yet ;)

Solution 5

Sobel gradient filter generator

(This answer refers to the analysis given by @Daniel, above.)

Gx[i,j] = i / (i*i + j*j)

Gy[i,j] = j / (i*i + j*j)

This is an important result, and a better explanation than can be found in the original paper. It should be written up in Wikipedia, or somewhere, because it also seems superior to any other discussion of the issue available on the internet.

However, it is not actually true that integer-valued representations are impractical for filters of size greater than 5*5, as claimed. Using 64-bit integers, Sobel filter sizes up to 15*15 can be exactly expressed.

Here are the first four; the result should be divided by the "weight", so that the gradient of an image region such as the following, is normalized to a value of 1.

1 2 3 4 5
1 2 3 4 5
1 2 3 4 5
1 2 3 4 5
1 2 3 4 5

Gx(3) :

-1/2  0/1  1/2           -1  0  1
-1/1    0  1/1   * 2 =   -2  0  2
-1/2  0/1  1/2           -1  0  1

weight = 4               weight = 8

Gx(5) :

-2/8 -1/5  0/4  1/5  2/8             -5  -4   0   4   5
-2/5 -1/2  0/1  1/2  2/5             -8 -10   0  10   8
-2/4 -1/1    0  1/1  2/4   * 20 =   -10 -20   0  20  10
-2/5 -1/2  0/1  1/2  2/5             -8 -10   0  10   8
-2/8 -1/5  0/4  1/5  2/8             -5  -4   0   4   5

weight = 12                          weight = 240

Gx(7) :

-3/18 -2/13 -1/10   0/9  1/10  2/13  3/18             -130 -120  -78    0   78  120  130
-3/13  -2/8  -1/5   0/4   1/5   2/8  3/13             -180 -195 -156    0  156  195  180
-3/10  -2/5  -1/2   0/1   1/2   2/5  3/10             -234 -312 -390    0  390  312  234
 -3/9  -2/4  -1/1     0   1/1   2/4   3/9   * 780 =   -260 -390 -780    0  780  390  260
-3/10  -2/5  -1/2   0/1   1/2   2/5  3/10             -234 -312 -390    0  390  312  234
-3/13  -2/8  -1/5   0/4   1/5   2/8  3/13             -180 -195 -156    0  156  195  180
-3/18 -2/13 -1/10   0/9  1/10  2/13  3/18             -130 -120  -78    0   78  120  130

weight = 24                                           weight = 18720

Gx(9) :

-4/32 -3/25 -2/20 -1/17  0/16  1/17  2/20  3/25  4/32                -16575  -15912  -13260   -7800       0    7800   13260   15912   16575
-4/25 -3/18 -2/13 -1/10   0/9  1/10  2/13  3/18  4/25                -21216  -22100  -20400  -13260       0   13260   20400   22100   21216
-4/20 -3/13  -2/8  -1/5   0/4   1/5   2/8  3/13  4/20                -26520  -30600  -33150  -26520       0   26520   33150   30600   26520
-4/17 -3/10  -2/5  -1/2   0/1   1/2   2/5  3/10  4/17                -31200  -39780  -53040  -66300       0   66300   53040   39780   31200
-4/16  -3/9  -2/4  -1/1     0   1/1   2/4   3/9  4/16   * 132600 =   -33150  -44200  -66300 -132600       0  132600   66300   44200   33150
-4/17 -3/10  -2/5  -1/2   0/1   1/2   2/5  3/10  4/17                -31200  -39780  -53040  -66300       0   66300   53040   39780   31200
-4/20 -3/13  -2/8  -1/5   0/4   1/5   2/8  3/13  4/20                -26520  -30600  -33150  -26520       0   26520   33150   30600   26520
-4/25 -3/18 -2/13 -1/10   0/9  1/10  2/13  3/18  4/25                -21216  -22100  -20400  -13260       0   13260   20400   22100   21216
-4/32 -3/25 -2/20 -1/17  0/16  1/17  2/20  3/25  4/32                -16575  -15912  -13260   -7800       0    7800   13260   15912   16575

weight = 40                                                          weight = 5304000

The Ruby program appended below, will calculate Sobel filters and corresponding weights of any size, although the integer-valued filters are not likely to be useful for sizes larger than 15*15.

#!/usr/bin/ruby

# Sobel image gradient filter generator
# by <[email protected]> -- Sept 2017
# reference:
# https://stackoverflow.com/questions/9567882/sobel-filter-kernel-of-large-size


if (s = ARGV[0].to_i) < 3 || (s % 2) == 0
    $stderr.puts "invalid size"
    exit false
end

s /= 2


n = 1

# find least-common-multiple of all fractional denominators
(0..s).each { |j|
    (1..s).each { |i|
        d = i*i + j*j
        n = n.lcm(d / d.gcd(i))
    }
}


fw1 = format("%d/%d", s, 2*s*s).size + 2
fw2 = format("%d", n).size + 2


weight = 0
s1 = ""
s2 = ""

(-s..s).each { |y|
    (-s..s).each { |x|
        i, j = x, y   # "i, j = y, x" for transpose
        d = i*i + j*j
        if (i != 0)
            if (n * i % d) != 0   # this should never happen
                $stderr.puts "inexact division: #{n} * #{i} / ((#{i})^2 + (#{j})^2)"
                exit false
            end
            w = n * i / d
            weight += i * w
        else
            w = 0
        end
        s1 += "%*s" % [fw1, d > 0 ? "%d/%d" % [i, d] : "0"]
        s2 += "%*d" % [fw2, w]
    }
    s1 += "\n" ; s2 += "\n"
}


f = n.gcd(weight)

puts s1

puts "\nweight = %d%s" % [weight/f, f < n ? "/%d" % (n/f) : ""]

puts "\n* #{n} =\n\n"

puts s2

puts "\nweight = #{weight}"

View more solutions

Comments

[Aarkan]

Thanks, any link to its generator function?

[Paul R]

I don't think these is one. Sobel is really only defined for 3x3 and the larger kernels seem to be determined on an ad hoc basis. They are just differentiators, so it should be fairly easy to generate coefficients for any size kernel.

[Jim G.]

But why is the scaling factor 1/8th for a 3x3 kernel and 1/16 for 5x5?

[Adam Bowen]

The 1D finite differences kernel is 1/2 * [1 0 -1] and the 1D smoothing kernel is 1/4 * [ 1 2 1 ]. So the 3x3 needs a scale of 1/2 * 1/4 = 1/8, and the larger kernels need 1/4 * 1/4 = 1/16 when moving from 1 to 2 dimensions.

[Daniel]

Sorry, but these are no Sobel kernels.

[Paul R]

@Daniel: would you care to expand on that ?

[Daniel]

Check my answer, it took a while to write down :-)

[Paul R]

Thanks for the comprehensive write-up - this is obviously a lot more rigorous than the approximations from the NIH document I based my answer on.

[Daniel]

If you're still interested in the topic, I would be grateful for feedback on my solution for Sobel kernels of arbitrary size and rotation.

[JasonN]

@Daniel I get a different answer for this. Following from Sobel's description, I'm unsure how all the roots fall out, and why you have different weightings for G(1,2)=-4 and G(2,1)=-8. They are both equidistant from the center where the derivative is being measured, and should have the same weight, no? Here's is what I got when I tried to derive NxN 2D Sobel operator, maybe I went wrong somewhere: reddit.com/r/computervision/comments/8e234p

[JasonN]

There's aren't even close.

[Paul R]

@Aarkan: please accept Daniel’s answer instead of this one.

[Daniel]

Hi! The -4 and -8 are not due to different distance weightings, but due to what I called i and j, the distance in gradient direction from the center element. If the entries were equal, you would measure diagonal gradients (in that case, a horizontal and vertical border would elicit the same response strength from the kernel). The second root factor should come from weighting pixel contributions with inverse distance, the first coming from normalizing the length of the "gradient contribution vector".

[JasonN]

Oh oh. " distance in gradient direction..." I see what you are saying. And hence why the center row/col is zeros (I want the projection onto the gradient direction). I was able to partially successfully rederive the 3x3 case by doing all the math by hand, but instead of 1's and 2's, I have sqrt(2)'s and 2's. Is that just for approximation or am I missing a step? When I tried to divide them out, everything cancels and I come up with Prewitt instead. Totally thanks. I was trying to make small changes on it for something else, and this really helped.

[Daniel]

I will expand the explanatory section for the 3x3 kernel, the details are too long for a comment

[Cris Luengo]

This is an interesting concept, since it follows Sobel's original derivation, but I'm not sure the resulting filter has good properties. You can make the kernel size very large, and see that it will not smooth much more than a small kernel (as entries further out approximate zero). Also, it is no longer separable, as the original Sobel kernel is.

[Daniel]

One thing to note though: If you have a uniform gradient in an image, the lower distance factor will just balance out the fact that the value difference is proportionally larger the farther out you go in gradient direction. In that special case, all pixels in gradient direction contribute the same amount. In the more common case of noisy or not-quite uniform gradients, larger kernels actually improve smoothing. For my purposes, 5x5 has so far usually been enough though.

[TylerH]

Hi @Daniel - would you mind voting on the synonymization request here? I'm trying to get [edge] synonymized with [edge-detection] as they are about the same thing, and users keep asking questions about microsoft-edge using the [edge] tag...

[eudoxos]

Thanks for inspiration, here are the routines (both integer and floating version) in c++: gist.github.com/eudoxos/b88f36e27eb782c6c60f25481c642aaf