1 Edge Detection

1.1 Introduction of Intensity Images

• An intensity image is a data matrix, whose values represent intensities within some range.

• Intensity image is represented as a single matrix, with each element of the matrix corresponding to one image pixel.

• The elements in the intensity matrix represent various intensities (or gray levels), where the intensity 0 usually represents 0 intensity (or black) and the intensity 1, 255, or 65535 usually represents full intensity (or white).

In matlab: To display an intensity image, use the imagesc (“image scale”) function.

1.2 Indexed Images

• An indexed image consists of a data matrix $X$, and a colormap matrix $map$.

• $map$ is an $m \times 3$ array of class double containing floatingpoint values in the range $[0, 1]$.

• Each row of $map$ specifies the red, green, and blue components of a single color.

1.3 Intensity Gradients

• The image is a function mapping coordinates to intensity $f(x,y)$

• The gradient of the intensity is a vector, which has an $x$ and a $y$ component. $$ \vec{G}[f(x,y)] = \left[ \begin{matrix} G_x \\ G_y \end{matrix} \right] = \left[ \begin{matrix} \frac{df}{dx} \\ \frac{df}{dy} \end{matrix} \right] $$

• Magnitude - $M(\vec{G}) = \sqrt{G^2_x + G^2_y}$

• Direction - $\alpha(x,y) = \mathbf{arctan}(\frac{G_y}{G_x})$

1.4 Approximating the Gradient

In image, it is diffcult to calculate the real gradient, but we can easily approximate the gradient.

For example: $$ G_x \cong f(i,j+1) - f(i,j) = \left[ \begin{matrix} -1 & 1 \end{matrix} \right] $$ $$ G_y \cong f(i,j) - f(i+1,j) = \left[ \begin{matrix} 1 \\ -1 \end{matrix} \right] $$

But this way calculates the gradient in different place. As shown below:

Thus, if We want to estimated in the same place, we can use a 2x2 mask instead. $$ G_x = \left[ \begin{matrix} -1 & 1 \\ -1 & 1 \end{matrix} \right] ,\ G_y = \left[ \begin{matrix} 1 & 1 \\ -1 & -1 \end{matrix} \right] $$

For each mask of weights, you multiply the corresponding pixel by the weight and sum over all pixels.

1.4.1 Other edge detectors

Roberts

Sobel

1.5 Convolution

Convolution is the computation of weighted sums of image pixels. For each pixel $p[i,j]$ in the image, the value $h[i,j]$ is calculated by translating the mask to pixel $p[i,j]$ and taking the weighted sum of pixels in neighbourhood.

1.6 Edge Detection

Steps:

• Take image
• Convolve mask with image for each direction ($G_x$ and $G_y$)
• Calculate magnitude $M(\vec{G}) = \sqrt{G^2_x + G^2_y}$. However, the cost of square root ($\sqrt{x}$) is expensive. Thus, it is normally replaced by $|G_x| + |G_y|$.

Note: We could detect edges by calculating the intensity change (gradient) across the image, and implement this using the idea of filtering.

2 Noise Filtering

We can use convolution to remove noise. Note: we should use noise filters before edge detectors.

2.1 Mean filter

One of filters is mean filter

Where, $n$ is the size of the filter. At here, the size of the filter is $3 \times 3$, so $n=9$. Obviously, mean filter is a linear filter.

% load data
shakey = read_image('','shakey.150.gif');
show_image(shakey);
title("Original")
load sobel.mat % for details, see Lab 1 on this blog

% Do not use nosie Noise filter, threshold is 30, using horizontal Sobel operator

shakey_X = conv2(shakey,sobelX,'valid');
show_image(abs(shakey_X)>30);
title("Gx")

% when size is 3 by 3, threshold is 30, using horizontal Sobel operator
sizeFilter = 3;
meanFilter = ones(sizeFilter)./(sizeFilter^2);

shakey_meanFilter = conv2(shakey,meanFilter,'valid');
show_image(shakey_meanFilter);
title("3by3")

shakey_X = conv2(shakey_meanFilter,sobelX,'valid');
show_image(abs(shakey_X)>30);
title("Gx 3by3")

% when size is 5 by 5, threshold is 30, using horizontal Sobel operator
sizeFilter = 5;
meanFilter = ones(sizeFilter)./(sizeFilter^2);

shakey_meanFilter = conv2(shakey,meanFilter,'valid');
show_image(shakey_meanFilter);
title("5by5")

shakey_X = conv2(shakey_meanFilter,sobelX,'valid');
show_image(abs(shakey_X)>30);
title("Gx 5by5")

2.2 Gaussian filter

Gaussian filtering is The most widely used.

In addition, We can replace a 2D Gaussian filter with two 1D Gaussian filters in sequence.

For example:
$$ \left[ \begin{matrix} 0.003 & 0.133 & 0.0219 & 0.133 & 0.003 \\ 0.133 & 0.0596 & 0.0983 & 0.0596 & 0.133 \\ 0.0219 & 0.0983 & 0.1621 & 0.0983 & 0.0219 \\ 0.133 & 0.0596 & 0.0983 & 0.0596 & 0.133 \\ 0.003 & 0.133 & 0.0219 & 0.133 & 0.003 \end{matrix} \right] $$

whcih can be decomposed into:
$$ \left[ \begin{matrix} 0.0545 & 0.2442 & 0.4026 & 0.2442 & 0.0545 \end{matrix} \right] \times \left[ \begin{matrix} 0.0545 \\ 0.2442 \\ 0.4026 \\ 0.2442 \\ 0.0545 \end{matrix} \right] $$

3 Advanced Edge Detection

3.1 Laplacian Operator (Second order operators)

In two dimensions, the definition of Laplacian Operator is: $$ \nabla^2 f = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} $$

As mentioned before, We approximate the derivative of $f$ in the $x$ direction as: $$ G_x = \frac{\partial f}{\partial x} = f(i,j+1) - f(i,j) $$

where, $f(i,j)$ corresponds to the coordinates of the pixel.

Thus, we have: $$ \begin{aligned} \frac{\partial^2 f}{\partial x^2} &= \frac{\partial G_x}{\partial x} \\ &= \frac{\partial [f(i,j+1) - f(i,j)]}{\partial x} \\ &= \frac{\partial f(i,j+1)}{\partial x} - \frac{\partial f(i,j)}{\partial x} \\ &= [f(i,j+2) - f(i,j+1)] - [f(i,j+1) - f(i,j)] \\ &= f(i,j+2) - 2*f(i,j+1) + f(i,j) \end{aligned} $$

Similarly, we can get $$ \frac{\partial^2 f}{\partial y^2} = f(i+2,j) - 2*f(i+1,j) + f(i,j) $$

Then, we get the Laplacian Operator as: $$ \nabla ^2 = \left[ \begin{matrix} 2 & -2 & 1 \\ -2 & 0 & 0 \\ 1 & 0 & 0 \end{matrix} \right] $$

However, in practice, the Laplacian operator is normally approximate as: $$ \nabla ^2 \approx \left[ \begin{matrix} 0 & -1 & 0 \\ -1 & 4 & -1 \\ 0 & -1 & 0 \end{matrix} \right] $$

Another commonly used small kernels are $$ \nabla ^2 \approx \left[ \begin{matrix} -1 & -1 & -1 \\ -1 & 8 & -1 \\ -1 & -1 & -1 \end{matrix} \right] $$

3.2 Laplacian of Gaussian (LoG)

We know the Laplacian Operator is: $$ L(x,y) = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} $$

The 2D Gaussian Distribution (mean = 0) function is: $$ G_\sigma(x,y) = \frac{1}{2 \pi \sigma^2} \cdot e^{-\frac{x^2+y^2}{2 \sigma^2}} $$

PS: The 1D Gaussian Distribution (mean = 0) function is: $$ G_\sigma(x) = \frac{1}{\sqrt{2 \pi} \sigma} \cdot e^{-\frac{x^2}{2 \sigma^2}} $$

When mean is 0, standard deviation is $\sigma$, we have: $$ LoG(x,y) = -\frac{1}{\pi \sigma^4} \cdot (1 - \frac{x^2+y^2}{2 \sigma^2}) \cdot e^{-\frac{x^2+y^2}{2 \sigma^2}} $$

Formula Derivation:

First, we know $$ LoG(x,y) = \frac{\partial^2 G_\sigma(x,y)}{\partial x^2} + \frac{\partial^2 G_\sigma(x,y)}{\partial y^2} $$

then, we solve $\frac{\partial^2 G_\sigma(x,y)}{\partial x^2}$ and $\frac{\partial^2 G_\sigma(x,y)}{\partial y^2}$ separately.

and

Similarly, we can get $$ \frac{\partial^2}{\partial y^2}G_\sigma(x,y) = \frac{1}{2 \pi \sigma^2} \cdot \frac{y^2 - \sigma^2}{\sigma^4} \cdot e^{-\frac{x^2+y^2}{2 \sigma^2}} $$

Thus, we get $$ \begin{aligned} LoG(x,y) &= \frac{\partial^2 G_\sigma(x,y)}{\partial x^2} + \frac{\partial^2 G_\sigma(x,y)}{\partial y^2} \\ &= \frac{1}{2 \pi \sigma^2} \cdot \frac{x^2 + y^2 -2 \sigma^2}{\sigma^4} \cdot e^{-\frac{x^2+y^2}{2 \sigma^2}} \\ &= -\frac{1}{\pi \sigma^4} \cdot (1 - \frac{x^2+y^2}{2 \sigma^2}) \cdot e^{-\frac{x^2+y^2}{2 \sigma^2}} \end{aligned} $$

For example:

% load data
shakey = read_image('','shakey.150.gif');

% creat filter
size = 9;
valueRange = -4:1:4;
s = 1.4; % standard deviation

LoG_Filter = zeros(size);

% calculate the value of LoG filter
for i = 1:size
    x = valueRange(i);
    LoG_Filter(i,:) = LoG(s,x,valueRange);
end

% apply LoG filter
shakeyLoG = conv2(shakey,LoG_Filter,'valid');

figure();
show_image(shakeyLoG);
title("shakeyLoG")

%% LoGfunction

function VLoG = LoG(s,x,y)

VLoG = -1/(s^4*pi).*(1 - (x.^2+y.^2)./(s^2*2)) .* exp(-(x.^2+y.^2)./(s^2*2));

end

Output:

3.3 Canny Edge Detector

Canny has shown that the first derivative of the Gaussian closely approximates the operator that optimizes the product of signal-to-noise ratio and localization. (analysis based on “step-edges” corrupted by “Gaussian noise”)

Steps of Algorithm:

1. Compute $f_x$ and $f_y$:

where $G(x,y)$ is the Gaussian function, $G_x(x,y)$ is the derivate of $G(x,y)$ with respect to $x$: $$ \begin{aligned} G_x(x,y) &= \frac{\partial}{\partial x} (\frac{1}{2 \pi \sigma^2} \cdot e^{-\frac{x^2+y^2}{2 \sigma^2}}) \\ &= \frac{1}{2 \pi \sigma^2} \cdot \frac{\partial}{\partial x}e^{-\frac{x^2+y^2}{2 \sigma^2}} \\ &= - \frac{1}{2 \pi \sigma^2} \cdot \frac{x}{\sigma^2} \cdot e^{-\frac{x^2+y^2}{2 \sigma^2}} \\ &= - \frac{x}{\sigma^2} \cdot G(x,y) \end{aligned} $$

in the same way, $G_y(x,y)$ is the derivate of $G(x,y)$ with respect to $y$: $$ G_y(x,y) = - \frac{y}{\sigma^2} \cdot G(x,y) $$

$\sigma$ is the standard deviation of $G(x,y)$.

2. Compute the gradient magnitude and direction: $$ magn(x,y) = |f_x| + |f_y| \\ dir(x,y) = \arctan(f_y/f_x) $$

3. Apply non-maxima suppression.

4. Apply hysteresis thresholding / edge linking.

4 Hough Transform

4.1 Hough Transform Basics

Assuming there is a point $p(x,y)$ in Cartesian Space, so Hough Transform provides a method to map $p(x,y)$ from Cartesian Space to Polar Space.

In Polar Space, we use $w$ and $\theta$ to define points, where $w$ is the distance and $\theta$ is angle from the origin (or reference point). Thus, we map $p(x,y)$ in Cartesian Space to $p(w,\theta)$ in Polar Space. $$ w = \sqrt{x^2 + y^2} \\ \theta = \arctan(y/x) $$

However, there are infinite lines can go through $p(x,y)$. How to define these lines?

Setting a reference point $p_o$, a point $p(x,y)$ located on a line $y = kx + b$. We use $w$ to express the distance between $p_o$ and this line, $\theta$ to express the angle between axis and the normal of this line.

Then, we use $-\frac{\cos\theta}{\sin\theta}$ to express $k$, $\frac{w}{\sin\theta}$ to express $b$. So that we have $$ \begin{aligned} y &= -\frac{\cos\theta}{\sin\theta} \cdot x + \frac{w}{\sin\theta} \\ y \sin\theta &= - x \cos\theta + w \\ w &= x \cos\theta + y \sin\theta \end{aligned} $$

Now, For a given point $p_1(x,y)$, we can calculate all lines going through it by function $w = x \cos\theta + y \sin\theta$.

4.2 Application

Normally, we use Hough Transform in Feature Extraction (特征提取). Hough Transform is used to identify features in objects, such as lines.

1. Each straight line in the image space corresponds to a single point in the Polar space;
2. Any part of the line segment on the straight line in the image space corresponds to the same point in the Polar space.

Thus, the line detection task is accomplished by finding peaks in the Polar space.

Lab 1

Default setting

1 Default functions

read_image.m

function image = read_image(base_dir, file)

% read_image(base_dir, file) - reads a gif file from the directory base_dir
%                        base_dir should be a string, and if it is the
%                        current directory it should be ''
% It also converts the image from unit8 to double

image = imread([base_dir, file]);

image = double(image);

show_image.m

function show_image(img)
% SHOW_IMAGE Image

figure % creates a new Figure window
colormap(gray);
imagesc(img);

2 Default parameters

sobel detectors

Roberts detectors

Tasks

Task 1: Create several different thresholds

Question:

What do you notice regarding the effect of changing the threshold?

Source Code

load sobel.mat

m = magnitude(sobelX,sobelY);

show_image(m);
title("m")

for i =20:10:40
show_image(m>i);
title("thresholds=" + i)
end


function m = magnitude(Gx,Gy)

shakey = read_image('','shakey.150.gif');

shakey_X = conv2(shakey,Gx,'valid');
shakey_Y = conv2(shakey,Gy,'valid');

m = (shakey_X.^2 + shakey_Y.^2).^(0.5);

end

Output

Task 2: Load up the Roberts operator

Question:

What do you notice regarding the difference between Sobel and Roberts

Source Code

load roberts.mat

m = magnitude(robertsA,robertsB);

show_image(m);
title("m")

for i =20:10:40
show_image(m>i);
title("thresholds=" + i)
end


function m = magnitude(Gx,Gy)

shakey = read_image('','shakey.150.gif');

shakey_X = conv2(shakey,Gx,'valid');
shakey_Y = conv2(shakey,Gy,'valid');

m = (shakey_X.^2 + shakey_Y.^2).^(0.5);

end

Output

Task 3: Replace the gradients with the absolute value of |Gx|+|Gy|

Question:

What do you notice regarding the difference between magnitude and absolute when calculating the edge?

Source Code

• Note: At here, I use roberts detectors.

load roberts.mat

m = magnitude(robertsA,robertsB);

show_image(m);
title("m")

for i =20:10:40
show_image(m>i);
title("thresholds=" + i)
end


function m = magnitude(Gx,Gy)

shakey = read_image('','shakey.150.gif');

shakey_X = conv2(shakey,Gx,'valid');
shakey_Y = conv2(shakey,Gy,'valid');

m = abs(shakey_X) + abs(shakey_Y);

end

Output

Lab 2

Default setting

1 Default functions

read_image.m

function image = read_image(base_dir, file)

% read_image(base_dir, file) - reads a gif file from the directory base_dir
%                        base_dir should be a string, and if it is the
%                        current directory it should be ''

image = imread([base_dir, file]);

image = double(image);

show_image.m

function show_image(img)
% SHOW_IMAGE Image

colormap(gray);
imagesc(img);

N.m

function p = N(m,s,x)

% N(m,s,x) - gives the Probability Density Function (pdf) of the Normal 
% Distribution with mean m, Standard Deviation (sd) s, for the value x.
% example: p = N(0,0.1,[-3:1:3]);
% Calculates pdf with mean = 0, sd = 0.1, for a vector of 7 elements long

p = 1/(s * sqrt(2* pi)) * exp(-(x-m).^2/(2*s^2));

2 Default parameters

sobel detectors

Tasks

Task 1

Question:

Can you describe the effect in comparison with applying the edge filter to the image directly?

Source Code

% load data
shakey = read_image('','shakey.150.gif');
load filters.mat


% The edge image without Gaussian noise filter, thresholds = 30
shakeyM = magnitude(shakey,sobelX,sobelY);
figure();
show_image(abs(shakeyM)>30);
title("Original edge image")

% The edge image with 3 by 3 Gaussian noise filter, thresholds = 30
shakeyGaussian3 = conv2(shakey,gaussian_filter_3x3,'valid');
shakeyM3 = magnitude(shakeyGaussian3,sobelX,sobelY);
figure();
show_image(abs(shakeyM3)>30);
title("gaussian filter 3x3")

% The edge image with 5 by 5 Gaussian noise filter, thresholds = 30
shakeyGaussian5 = conv2(shakey,gaussian_filter_5x5,'valid');
shakeyM5 = magnitude(shakeyGaussian5,sobelX,sobelY);
figure();
show_image(abs(shakeyM5)>30);
title("gaussian filter 5x5")


% Calculating edge
function m = magnitude(img,Gx,Gy)

shakey_X = conv2(img,Gx,'valid');
shakey_Y = conv2(img,Gy,'valid');

m = abs(shakey_X) + abs(shakey_Y);
end

Output

Where I found that noise level dropped.

Task 2

Question:

What is the effect of increasing the size of the Gaussian Filter (3x3 versus 5x5 for example)? What is the effect of changing the standard deviation s? Why do you see what you see?

Source Code

% load data
shakey = read_image('','shakey.150.gif');
load filters.mat

% The edge image without Gaussian noise filter, thresholds = 30
shakeyM = magnitude(shakey,sobelX,sobelY);
figure();
show_image(abs(shakeyM)>30);
title("Original edge image")

% Creat gaussian filter
m = 0; % mean
s = 5; % standard deviation
ValueCoordinates = -3 : 1.5 : 3;
GaussianCreated = N(m,s,ValueCoordinates);
GaussianFilter5 = GaussianCreated'*GaussianCreated;

ValueCoordinates = -6 : 1.5 : 6;
GaussianCreated = N(m,s,ValueCoordinates);
GaussianFilter9 = GaussianCreated'*GaussianCreated;

% The edge image with 5 by 5 Gaussian noise filter, thresholds = 30
shakeyGaussian5 = conv2(shakey,GaussianFilter5,'valid');
shakeyM5 = magnitude(shakeyGaussian5,sobelX,sobelY);

% The edge image with 9 by 9 Gaussian noise filter, thresholds = 30
shakeyGaussian9 = conv2(shakey,GaussianFilter9,'valid');
shakeyM9 = magnitude(shakeyGaussian9,sobelX,sobelY);

figure();
show_image(abs(shakeyM5)>30);
title("gaussian filter 5x5, edge")
figure();
show_image(shakeyGaussian5);
title("gaussian filter 5x5, standard deviation is " + s)

figure();
show_image(abs(shakeyM9)>30);
title("gaussian filter 9x9, edge")
figure();
show_image(shakeyGaussian9);
title("gaussian filter 9x9, standard deviation is " + s)



%% Calculating edge
function m = magnitude(img,Gx,Gy)

shakey_X = conv2(img,Gx,'valid');
shakey_Y = conv2(img,Gy,'valid');

m = abs(shakey_X) + abs(shakey_Y);

end

Output

Where I found that with the increasing of standard deviation, the images get blurrier and the gradient between pixels decreases.

According to the above pictures, I found the size of the Gaussian Filter should be changed with the change of standard deviation. The big value of standard deviation corresponds to big size of the Gaussian Filter.

Task 3

Question:

I mentioned the Laplacian of the Gaussian in the lecture. How could you combine the idea of the Laplacian operator with the idea of Gaussian smoothing? Try out your ideas.

Source Code

% load data
shakey = read_image('','shakey.150.gif');

% creat filter
size = 9;
valueRange = -4:1:4;
s = 1.4; % standard deviation

LoG_Filter = zeros(size);

% calculate the value of LoG filter
for i = 1:size
    x = valueRange(i);
    LoG_Filter(i,:) = LoG(s,x,valueRange);
end

% apply LoG filter
shakeyLoG = conv2(shakey,LoG_Filter,'valid');

figure();
show_image(shakeyLoG);
title("shakeyLoG")

%% LoGfunction

function VLoG = LoG(s,x,y)

VLoG = -1/(s^4*pi).*(1 - (x.^2+y.^2)./(s^2*2)) .* exp(-(x.^2+y.^2)./(s^2*2));

end

Output

Source Code

load data
shakey = read_image('','shakey.150.gif');
load filters.mat

% Creat gaussian filter
m = 0; % mean
s = 1; % standard deviation
ValueCoordinates = -3 : 1 : 3;
GaussianCreated = N(m,s,ValueCoordinates);
GaussianFilter7 = GaussianCreated'*GaussianCreated;

% The edge image with 7 by 7 Gaussian noise filter
shakeyGaussian7 = conv2(shakey,GaussianFilter7,'valid');

% use Sobel filter determine the edge
shakeyGaussian7Sobel = magnitude(shakeyGaussian7,sobelX,sobelY);

% use Laplacian filter determine the edge
shakeyGaussian7Laplacian = edge(shakeyGaussian7,'zerocross');

figure();
show_image(abs(shakeyGaussian7Sobel)>30);
title("Gaussian7Sobel, standard deviation is 1, thresholds = 30")

figure();
show_image(shakeyGaussian7Laplacian);
title("Gaussian7Laplacian, standard deviation is 1")



%% Calculating edge
function m = magnitude(img,Gx,Gy)

shakey_X = conv2(img,Gx,'valid');
shakey_Y = conv2(img,Gy,'valid');

m = abs(shakey_X) + abs(shakey_Y);

end

Output

We can find that the outputs are totally different between LoG and the way that using Gaussian and Laplacianseparately.