Image Processing - Implementing Sobel Filter

Question

I've got a task to implement Sobel filter which is, as you know, an image processing filter for edge detection. But unfortunately, I've got no experience in image processing field, to the extent that I don't even know how images are represented in computer. Totally no knowledge in this field.

I've read some papers and PDFs but they focus on many topics which I feel that I may not need them for my task.

I'd be happy to know your suggestions or if there is any particular paper, PDF, tutorial or quick guide for this purpose.

Thank you

EDIT:

Thank you all :) The result of our work can be downloaded from here.

Answer 1

It's pretty easy, you just need to convolve your image with a Sobel filter. A Sobel filter has two kernels, x-direction kernel and y-direction kernel. The x-direction kernel detects horizontal edges, and y-direction kernels detects vertical edges.

x-direction kernel (the size is 3x3)

float kernelx[3][3] = {{-1, 0, 1}, 
                       {-2, 0, 2}, 
                       {-1, 0, 1}};

y-direction kernel

float kernely[3][3] = {{-1, -2, -1}, 
                        {0,  0,  0}, 
                        {1,  2,  1}};

To calculate the convolution at pixel (x,y), define a window of size equal to the kernel size (source code to calculate magnitude in x and magnitude in y are identical):

double magX = 0.0; // this is your magnitude

for(int a = 0; a < 3; a++)
{
    for(int b = 0; b < 3; b++)
    {            
        int xn = x + a - 1;
        int yn = y + b - 1;

        int index = xn + yn * width;
        magX += image[index] * kernelx[a][b];
    }
 }

Note that the input is a grayscale image and it can be represented as 1D array of double (This is just a trick, since a pixel value in coordinate (x,y) can be accessed with index = [x + y * width] )

To calculate magnitude in pixel (x,y) given magX and magY :

mag = sqrt( magX^2 + magY^2 )

Answer 2

The most simple explanation of the Sobel operator I've seen to this date is from Saush's blog, a tech enthusiast who once met Sobel himself:

The post describes in (not too many) details how to implement the filter, and shares Ruby source-code for demonstration purposes:

require 'chunky_png'

class ChunkyPNG::Image
  def at(x,y)
    ChunkyPNG::Color.to_grayscale_bytes(self[x,y]).first
  end
end

img = ChunkyPNG::Image.from_file('engine.png')

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

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

edge = ChunkyPNG::Image.new(img.width, img.height, ChunkyPNG::Color::TRANSPARENT)

for x in 1..img.width-2
  for y in 1..img.height-2
    pixel_x = (sobel_x[0][0] * img.at(x-1,y-1)) + (sobel_x[0][1] * img.at(x,y-1)) + (sobel_x[0][2] * img.at(x+1,y-1)) +
              (sobel_x[1][0] * img.at(x-1,y))   + (sobel_x[1][1] * img.at(x,y))   + (sobel_x[1][2] * img.at(x+1,y)) +
              (sobel_x[2][0] * img.at(x-1,y+1)) + (sobel_x[2][1] * img.at(x,y+1)) + (sobel_x[2][2] * img.at(x+1,y+1))

    pixel_y = (sobel_y[0][0] * img.at(x-1,y-1)) + (sobel_y[0][1] * img.at(x,y-1)) + (sobel_y[0][2] * img.at(x+1,y-1)) +
              (sobel_y[1][0] * img.at(x-1,y))   + (sobel_y[1][1] * img.at(x,y))   + (sobel_y[1][2] * img.at(x+1,y)) +
              (sobel_y[2][0] * img.at(x-1,y+1)) + (sobel_y[2][1] * img.at(x,y+1)) + (sobel_y[2][2] * img.at(x+1,y+1))

    val = Math.sqrt((pixel_x * pixel_x) + (pixel_y * pixel_y)).ceil
    edge[x,y] = ChunkyPNG::Color.grayscale(val)
  end
end

edge.save('engine_edge.png')

Input/Output:

Answer 3

Sobel Operator Wikipedia page is well descriptive about how to perform it. There other operators such as Roberts cross and Prewitt

Using convolution operation, you can switch the approach by changing the kernel matrix. Below, the implementation of Sobel and Convolution using Marvin Framework may help you.

Sobel:

public class Sobel extends MarvinAbstractImagePlugin{

    // Definitions
    double[][] matrixSobelX = new double[][]{
            {1,     0,  -1},
            {2,     0,  -2},
            {1,     0,  -1}
    };
    double[][] matrixSobelY = new double[][]{
            {-1,    -2,     -1},
            {0,     0,      0},
            {1,     2,      1}
    };

    private MarvinImagePlugin   convolution;

    public void load(){
        convolution = MarvinPluginLoader.loadImagePlugin("org.marvinproject.image.convolution.jar");
    }

    public MarvinAttributesPanel getAttributesPanel(){
        return null;
    }
    public void process
    (
        MarvinImage imageIn, 
        MarvinImage imageOut,
        MarvinAttributes attrOut,
        MarvinImageMask mask, 
        boolean previewMode
    )
    {
        convolution.setAttribute("matrix", matrixSobelX);
        convolution.process(imageIn, imageOut, null, mask, previewMode);
        convolution.setAttribute("matrix", matrixSobelY);
        convolution.process(imageIn, imageOut, null, mask, previewMode);
    }
}

Convolution:

public class Convolution extends MarvinAbstractImagePlugin{

    private MarvinAttributesPanel   attributesPanel;
    private MarvinAttributes        attributes;

    public void process
    (
        MarvinImage imageIn, 
        MarvinImage imageOut,
        MarvinAttributes attributesOut,
        MarvinImageMask mask, 
        boolean previewMode
    )
    {
        double[][] matrix = (double[][])attributes.get("matrix");

        if(matrix != null && matrix.length > 0){
            for(int y=0; y<imageIn.getHeight(); y++){
                for(int x=0; x<imageIn.getWidth(); x++){
                    applyMatrix(x, y, matrix, imageIn, imageOut);
                }
            }
        }
    }

    private void applyMatrix
    (
        int x,
        int y,
        double[][] matrix,
        MarvinImage imageIn,
        MarvinImage imageOut
    ){

        int nx,ny;
        double resultRed=0;
        double resultGreen=0;
        double resultBlue=0;

        int xC=matrix[0].length/2;
        int yC=matrix.length/2;

        for(int i=0; i<matrix.length; i++){
            for(int j=0; j<matrix[0].length; j++){
                if(matrix[i][j] != 0){      
                    nx = x + (j-xC);
                    ny = y + (i-yC);

                    if(nx >= 0 && nx < imageOut.getWidth() && ny >= 0 && ny < imageOut.getHeight()){

                        resultRed   +=  (matrix[i][j]*(imageIn.getIntComponent0(nx, ny)));
                        resultGreen +=  (matrix[i][j]*(imageIn.getIntComponent1(nx, ny)));
                        resultBlue  +=  (matrix[i][j]*(imageIn.getIntComponent2(nx, ny)));
                    }


                }



            }
        }

        resultRed   = Math.abs(resultRed);
        resultGreen = Math.abs(resultGreen);
        resultBlue = Math.abs(resultBlue);

        // allow the combination of multiple appications
        resultRed   += imageOut.getIntComponent0(x,y);
        resultGreen += imageOut.getIntComponent1(x,y);
        resultBlue  += imageOut.getIntComponent2(x,y);

        resultRed   = Math.min(resultRed, 255);
        resultGreen = Math.min(resultGreen, 255);
        resultBlue  = Math.min(resultBlue, 255);

        resultRed   = Math.max(resultRed, 0);
        resultGreen = Math.max(resultGreen, 0);
        resultBlue  = Math.max(resultBlue, 0);

        imageOut.setIntColor(x, y, imageIn.getAlphaComponent(x, y), (int)resultRed, (int)resultGreen, (int)resultBlue);
    }

    public void load(){
        attributes = getAttributes();
        attributes.set("matrix", null);
    }

    public MarvinAttributesPanel getAttributesPanel(){
        if(attributesPanel == null){
            attributesPanel = new MarvinAttributesPanel();
            attributesPanel.addMatrixPanel("matrixPanel", "matrix", attributes, 3, 3);
        }
        return attributesPanel;
    }

}

Answer 4

Gx is estimating the gradient in the x-direction (columns) and Gy is estimating the gradient in the y-direction (rows). So Gy detects horizontal lines, and Gx detects vertical lines.

Answer 5

Of course, you could use OpenCV for this:

import cv2
import numpy as np

img = cv2.imread(INPUT_IMAGE)
img = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY).astype(float)

edge_x = cv2.Sobel(img,cv2.CV_64F,1,0,ksize=3)
edge_y = cv2.Sobel(img,cv2.CV_64F,0,1,ksize=3)    
edge = np.sqrt(edge_x**2 + edge_y**2)    # image can be normalized to 
                                         # fit into 0..255 color space
cv2.imwrite(OUTPUT_IMAGE, edge)

Input / Output:

Answer 6

All the above mentioned steps in a R markdown file. Hopefully this makes it more visual and easier to understand. I needed to implement a sobel filter and this page helped me understand the concepts still I had some trouble in getting it done. So putting it all in one place hopefully it helps.

http://rpubs.com/ghub_24/420754

Answer 7

I use Octave 4.4.1 for image processing and edge detection. Octave is an open source software which provides functionality of MatLab. I have used the following code from "https://in.mathworks.com/" Let us consider that the image is 'k1.jpg' To implement Sobel edge operator i=imread('k1.jpg');//To read the image pkg load image//Loading image processing Tool Box g=rgb2gray(i)//Conversion of image to grayscale S=edge(g,'Sobel');//Applying Sobel edge detector to g imshow(S);//To display the image Original image Segmented Image

Answer 8

You can do this in R with the raster package (oriented towards geographic data)

library(raster)
sobel <- function(r) {
    fy <- matrix(c(1,0,-1,2,0,-2,1,0,-1)/4, nrow=3)
    fx <- matrix(c(-1,-2,-1,0,0,0,1,2,1)/4 , nrow=3)
    rx <- focal(r, fx)
    ry <- focal(r, fy)
    sqrt(rx^2 + ry^2)
}

b <- brick("https://i.stack.imgur.com/Bnxa6.jpg")
plotRGB(b)
s <- stack(lapply(1:3, function(i) sobel(b[[i]])))
plotRGB(s)

# or
ss <- mean(s)
plot(ss, col=gray(1:10/10))

# or
bb <- mean(b)
sss <- sobel(bb)
plot(sss, col=gray(1:10/10))