2

I'm trying to implement the sobel operator in Python and visualize it. However, I'm struggling with how to do that. I have the following code, which currently calculates the gradient for each pixel.

from PIL import Image
import math


def run():

    try:

        image = Image.open("brick-wall-color.jpg")
        image = image.convert('LA')

        apply_sobel_masks(image)

    except RuntimeError, e:
        print e


def apply_sobel_masks(image):

    gx = [
        [-1, 0, 1],
        [-2, 0, 2],
        [-1, 0, 1]
    ]

    gy = [
        [1, 2, 1],
        [0, 0, 0],
        [-1, -2, -1]
    ]

    width, height = image.size

    for y in range(0, height):

        for x in range(0, width):

            gradient_y = (
                gy[0][0] * get_pixel_safe(image, x - 1, y - 1, 0) +
                gy[0][1] * get_pixel_safe(image, x, y - 1, 0) +
                gy[0][2] * get_pixel_safe(image, x + 1, y - 1, 0) +
                gy[2][0] * get_pixel_safe(image, x - 1, y + 1, 0) +
                gy[2][1] * get_pixel_safe(image, x, y + 1, 0) +
                gy[2][2] * get_pixel_safe(image, x + 1, y + 1, 0)
            )

            gradient_x = (
                gx[0][0] * get_pixel_safe(image, x - 1, y - 1, 0) +
                gx[0][2] * get_pixel_safe(image, x + 1, y - 1, 0) +
                gx[1][0] * get_pixel_safe(image, x - 1, y, 0) +
                gx[1][2] * get_pixel_safe(image, x + 1, y, 0) +
                gx[2][0] * get_pixel_safe(image, x - 1, y - 1, 0) +
                gx[2][2] * get_pixel_safe(image, x + 1, y + 1, 0)
            )

            print "Gradient X: " + str(gradient_x) + " Gradient Y: " + str(gradient_y)
            gradient_magnitude = math.sqrt(pow(gradient_x, 2) + pow(gradient_y, 2))

            image.putpixel((x, y), #tbd)


    image.show()


def get_pixel_safe(image, x, y, layer):

    try:
        return image.getpixel((x, y))[layer]

    except IndexError, e:
        return 0


run()

Now the gradient_magnitude is often a value well outside of the 0-255 range e.g. 990.0, 1002.0, 778, etc.

So what I'd like to do is visualize that gradient, but I'm not sure how. Most resources online only mention calculating the gradient angle and magnitude but not how to represent it in an image.

saurabheights
  • 3,967
  • 2
  • 31
  • 50
fbailey
  • 303
  • 5
  • 17
  • 1
    Try normalizing to 0-255 value using the maximum and minimum value of gradient magnitude image, i.e. `displayableGradMagnitude = 255* (gradMagnitude-minGradMagnitude)/(maxGradMagnitude-minGradMagnitude)`. – saurabheights Mar 13 '17 at 21:16
  • Thank you, that worked perfectly. Maybe post that as the answer so I can give you your credit – fbailey Mar 13 '17 at 21:34
  • Also note that the Solbel convolution kernel is separable. You can apply a convolution with [-1,0,1], and on the result a second convolution with the transpose of [1,2,1]. That is always more efficient, and it's also easier to implement. The order of the two operations is irrelevant. The y-derivative is computed by transposing the two kernels. – Cris Luengo Mar 14 '17 at 00:47
  • @CrisLuengo, could you elaborate more? Say I have a 3x3 grid of pixels, I'm just multiplying pixel[0] by -1, pixel[1] by 0, pixel[3] by 1 and repeating for the following two rows? and then repeating once more with the second convolution? – fbailey Mar 14 '17 at 03:53
  • 2
    @fbailey: For each row, do `tmp[x]=in[x-1]-in[x+1]`. Then, for each column, do `out[y]=tmp[y-1]+2*tmp[y]+tmp[y+1]`. That gives you the x-derivative. Sorry for the speudo-code here, these comment boxes don't really let you do much more. But I'm sure you can figure it out from there. Note that you do not need to create a temporary image `tmp`, you can write directly in the output image for each of the two steps. But in the second step you then need to copy the column you'll be working on so you don't overwrite data you still need. There's cleverer ways of doing that, but that is a longer story. – Cris Luengo Mar 14 '17 at 04:02
  • 1
    @fbailey: You are using 3x3 pixels at each location to compute gradient, but what chris has mentioned is you can divide this processing into first X-direction and then Y-direction computation(dividing the 3x3 kernel into one horizontal and one vertical vector). Read DIP by Gonzalez, it has all this info. P.S. I will post answer within 12 hours. – saurabheights Mar 14 '17 at 04:09

2 Answers2

3

Using @saurabheights suggestion I was able to visualize the gradient's magnitude. One mistake I also corrected, was that I was editing each pixel after calculating its gradient. This is incorrect because when the kernel shifts over one pixel it now uses the value of the just edited pixel. The corrected code is posted below:

from PIL import Image, ImageFilter
import math


def run():

    try:

        image = Image.open("geo.jpg")
        image = image.convert('LA')
        image = image.filter(ImageFilter.GaussianBlur(radius=1))
        apply_sobel_masks(image)

    except RuntimeError, e:
        print e


def apply_sobel_masks(image):

    gx = [
        [-1, 0, 1],
        [-2, 0, 2],
        [-1, 0, 1]
    ]

    gy = [
        [1, 2, 1],
        [0, 0, 0],
        [-1, -2, -1]
    ]

    width, height = image.size
    gradient_magnitudes = [[0 for x in range(width)] for y in range(height)]
    gradient_max = None
    gradient_min = None

    for y in range(0, height):

        for x in range(0, width):

            gradient_y = (
                gy[0][0] * get_pixel_safe(image, x - 1, y - 1, 0) +
                gy[0][1] * get_pixel_safe(image, x, y - 1, 0) +
                gy[0][2] * get_pixel_safe(image, x + 1, y - 1, 0) +
                gy[2][0] * get_pixel_safe(image, x - 1, y + 1, 0) +
                gy[2][1] * get_pixel_safe(image, x, y + 1, 0) +
                gy[2][2] * get_pixel_safe(image, x + 1, y + 1, 0)
            )

            gradient_x = (
                gx[0][0] * get_pixel_safe(image, x - 1, y - 1, 0) +
                gx[0][2] * get_pixel_safe(image, x + 1, y - 1, 0) +
                gx[1][0] * get_pixel_safe(image, x - 1, y, 0) +
                gx[1][2] * get_pixel_safe(image, x + 1, y, 0) +
                gx[2][0] * get_pixel_safe(image, x - 1, y - 1, 0) +
                gx[2][2] * get_pixel_safe(image, x + 1, y + 1, 0)
            )

            gradient_magnitude = math.ceil(math.sqrt(pow(gradient_x, 2) + pow(gradient_y, 2)))

            if gradient_max is None:
                gradient_max = gradient_magnitude
                gradient_min = gradient_magnitude

            if gradient_magnitude > gradient_max:
                gradient_max = gradient_magnitude

            if gradient_magnitude < gradient_min:
                gradient_min = gradient_magnitude

            gradient_magnitudes[y][x] = gradient_magnitude

    # Visualize the gradients
    for y in range(0, height):

        for x in range(0, width):

            gradient_magnitude = gradient_magnitudes[y][x]
            pixel_value = int(math.floor(255 * (gradient_magnitude - gradient_min) / (gradient_max - gradient_min)))

            image.putpixel((x, y), pixel_value)

    image.show()


def get_pixel_safe(image, x, y, layer):

    try:
        return image.getpixel((x, y))[layer]

    except IndexError, e:
        return 0


run()
fbailey
  • 303
  • 5
  • 17
2

The simplest way to bring values into a specific range is to normalize. For n values, find minimum and maximum of all these values. For range [a, b], normalize each value x as:-

x' = a + (b-a) * (x-min)/(max-min)

For OP's scenario, this equation for gradient magnitude will be:-

x' = 255 * (x-min)/(max-min)

saurabheights
  • 3,967
  • 2
  • 31
  • 50