## Description

In this lab you will implement a differential geometry based edge detector that works on multiple

scales, and apply the results from this operator to detect lines with the Hough transform.

The goal of this lab is for you to understand how differential geometry based edge detection and

the Hough transform works as well as gain practical experience in applying these operations to real

data and learn the characteristics of the methods, including how important the choice of scale is

for edge detection.

As prerequisites to this lab you should have read the course material on differential operators,

edge detection and Hough transform. You should also preferably have completed and presented the

results from Lab 1.

Reporting: When reporting this lab emphasis is placed on: For parts 1-4 the experimental results

and interpretations/explanation of these, are the most central. For parts 5-6 you should also have

created well functioning and structured Matlab functions, that respectively perform edge detection

and accumulation in the Hough domain.

You are recommended to use the command subplot to assemble multiple results into figures, thus

simplifying the interpretation of these. In the answer sheet write down summarizing conclusions

and compile results for the explicitly stated exercises and questions.

Advice: Read through all the lab instructions before you begin and prepare some principle solutions

before you start writing the code to Exercise 6. Read through the hints and practical suggestions

in section 6.1. Perform some simple experiments to analyze the behavior for different coordinate

systems used to access pixels in images and curves respectively.

## 1 Difference operators

Create two difference operators deltax and deltay that approximate the first order partial derivatives

in two orthogonal directions. You may choose to use either the simple difference operator, central

differences, Robert’s diagonal operator or the Sobel operator. Then load the (function based) image

few256 from the course library with

tools = few256;

and compute discrete derivation approximations using

dxtools = conv2(tools, deltax, ’valid’);

dytools = conv2(tools, deltay, ’valid’);

Show the results on screen.

Question 1: What do you expect the results to look like and why? Compare the size of dxtools

with the size of tools. Why are these sizes different?

## 2 Point–wise thresholding of gradient magnitudes

Based on the results above, compute an approximation of the gradient magnitude

gradmagntools = sqrt(dxtoolsconv .^2 + dytoolsconv .^2);

and show the results on screen. Compute the histogram of this image and use this to guess a threshold

that yields reasonably thin edges when applied to gradmagntools. Show the thresholded image with

showgrey((gradmagntools – threshold) > 0)

for different values of threshold. Perform the same study on the image godthem256, whose content

include image structures of higher complexity.

The work is facilitated considerably if you write a

Matlab procedure

function pixels = Lv(inpic, shape)

if (nargin < 2)

shape = ’same’;

end

Lx = filter2(dxmask, inpic, shape);

Ly = filter2(dymask, inpic, shape);

pixels = Lx.^2 + Ly.^2;

where dxmask and dymask are Matlab procedures that return suitable masks that approximate

the derivatives. Preferably, combine these operations with smoothing of the image using a Gaussian

convolution step prior to computing the gradient magnitude. (Apply the procedure you developed in

previous lab, alternatively use discgaussfft.)

Question 2: Is it easy to find a threshold that results in thin edges? Explain why or why not!

Question 3: Does smoothing the image help to find edges?

## 3 Differential geometry based edge detection: Theory

One way of extracting thin edges is by considering points for which the gradient magnitude reaches

local maxima in gradient direction. As it has been mentioned during the lectures, this edge definition

can be expressed in differential geometry terms as the second order derivative being equal to zero and

the third order derivative being negative. Introduce in each point a local coordinate system (u, v) such

that the v direction is parallel to the gradient direction and the u direction is perpendicular to this.

Then an edge can be defined as

Lvv = 0,

Lvvv < 0,

where Lvv and Lvvv denote the second and third order derivatives of the smoothened intensity function

L in the v direction. Typically, we get L from the original image function f by convolving it with a

Gaussian kernel g

L(·; t) = g(·; t) ∗ f

where

g(x, y; t) = 1

2πt

e

−(x

2+y

2

)/(2t)

.

To express the edge definition above in partial derivatives of the smoothened intensity function L

with respect to the Cartesian coordinate directions, we write the gradient vector at an arbitrary point

(x0, y0) as

∇L =

Lx

Ly

(x0,y0)

= |∇L|

cos ϕ

sin ϕ

(x0,y0)

In terms of

cos ϕ = q

Lx

L2

x + L2

y

sin ϕ =

Ly

q

L2

x + L2

y

the derivatives in the u and v directions can be written as

∂u = sin ϕ ∂x − cos ϕ ∂y, ∂v = cos ϕ ∂x + sin ϕ ∂y.

After expansion the explicit expressions of Lvv and Lvvv assume the following forms

Lvv =

L

2

xLxx + 2LxLyLxy + L

2

yLyy

L2

x + L2

y

,

Lvvv = =

L

3

xLxxx + 3L

2

xLyLxxy + 3LxL

2

yLxyy + L

3

yLyyy

(L2

x + L2

y

)

3/2

.

Since only the sign of the derivatives are of interest for the edge definition, we can avoid dividing by

the gradient magnitude Lv = |∇L| =

q

L2

x + L2

y and instead define edges the following way:

L˜

vv = L

2

vLvv = L

2

xLxx + 2LxLyLxy + L

2

yLyy = 0,

L˜

vvv = L

3

vLvvv = L

3

xLxxx + 3L

2

xLyLxxy + 3LxL

2

yLxyy + L

3

yLyyy < 0.

## 4 Computing differential geometry descriptors

Write two Matlab procedures Lvvtilde and Lvvvtilde that compute the differential expressions L˜

vv

and L˜

vvv. Preferably, design this function in the same way as the Matlab procedure Lv in exercise 2.

You need to define masks corresponding to the discrete derivative approximations of all partial

derivatives up to the order of three. The easiest way to do this is by using the central differences of

the first order derivatives

δx =

−

1

2

, 0,

1

2

δy =

−

1

2

0

1

2

and simple difference approximations of the second order derivatives

δxx = (1, −2, 1) δyy =

1

−2

1

.

Approximations of higher order derivatives can then be created by applying these in cascade.

The notation here is a form of operators and means that every difference operator should be

interpreted as linear filtering with a mask. The concatenation of two operators thus becomes a filtering

operation of their corresponding masks. In Matlab this concatenation of two masks is performed with

the functions filter2 or conv2, with the argument ’shape’ set to ’same’. For example,

δxxx = δx δxx, δxy = δx δy, δxxy = δxx δy.

Note! When you create masks for discrete derivative approximations and apply these to image

data the convention depend on the definitions used for the coordinate axes, as well as on whether you

use the function filter2 or conv2.

Think through the definitions and be consistent

When you create masks for these computational molecules, think of how the argument shape of

the Matlab functions filter2 and conv2 treats those image points for which some of the points

of the filter masks will be placed outside the available amount of image data. In order for the image

boundaries to be treated consistently for all derivative approximations, you should let all filter masks

have the same size (e.g. 5 × 5), even if some masks then will contain a large number of zeros.

In order to make sure that your filter masks have the expected effect, you should analyze their

performances on some reference data. The easiest way to do this in Matlab is to define some matrices

with x and y coordinates according to

[x y] = meshgrid(-5:5, -5:5)

and study the results of for example the operations

filter2(dxxxmask, x .^3, ’valid’)

filter2(dxxmask, x .^3, ’valid’)

filter2(dxxymask, x .^2 .* y, ’valid’)

When these difference approximations are applied to polynomials the following should hold

δx(x

n

) = nxn−1 + (lower order terms), (1)

δxn (x

n

) = n! (2)

δxn+k (x

n

) = 0 (3)

δxn (y

k

) = 0 (4)

Make sure that this is the case.

Experiments: When you feel convinced that all the subcomponents work, load the image

house = godthem256;

and compute the zero crossings of L˜

vv on a number of different scales by writing

contour(Lvvtilde(discgaussfft(house, scale ), ’same’), [0 0])

axis(’image’)

axis(’ij’)

where you appropriately try scale = 0.0001, 1.0, 4.0, 16.0 and 64.0.

Question 4: What can you observe? Provide explanation based on the generated images.

Study the sign of the third order derivative in the gradient direction by loading the image

tools = few256;

and show the result of

showgrey(Lvvvtilde(discgaussfft(tools, scale), ’same’) < 0)

for the same values of scale. What is the effect of the sign condition in this differential expression?

Question 5: Assemble the results of the experiment above into an illustrative collage with the

subplot command. Which are your observations and conclusions?

Question 6: How can you use the response from L˜

vv to detect edges, and how can you improve

the result by using L˜

vvv?

## 5 Extraction of edge segments

So far we have only considered the result of computing differential operators at different scales and

looked at their zero crossings and sign variations. In this exercise we will collect these components to

a Matlab function

function edgecurves = extractedge(inpic, scale, threshold, shape)

that computes the edges of an image at arbitrary scale and returns these lists of edge points. If

the parameter threshold is given, this value shall be used as threshold for the gradient magnitude

computed at the same scale. Otherwise, no thresholding shall be used and all edge segments will be

returned.

Two functions are provided in the course library. The first one extracts level curves in a given image

and rejects points based on the sign of the second input argument. The second function thresholds

these curves with respect to the sign of another image.

function curves = zerocrosscurves(zeropic, maskpic1)

function curves = thresholdcurves(curves, maskpic2)

These functions are based on the Matlab level curve algorithm and Matlab’s internal format for

representation of polygons (see help contourc regarding the curve format, and respectively help

zerocrosscurves and help thresholdcurves for the functions of the subroutines in question). Even

if this format does not permit the best possible abstraction level, it has the advantage that it can

readily be used by the functions already existing in Matlab.

Write the function extractedge !

The exercise thus consists of, based on the theory described in section 3 and the analogy with the

experiments you performed in section 4, computing suitable differential quantities that can be used

as input data for the functions zerocrosscurves and thresholdcurves, so as to determine edges

thresholded with respect to the gradient magnitude at arbitrary scale. When your function works

satisfactory, use it to extract edges from the images house and tools. First try different scales e.g.

0.0001, 1.0, 4.0, 16.0 and 64.0. For each image find the best scale and adjust the threshold accordingly.

Preferably, show the results on screen with

overlaycurves(image, curves)

that overlays the curves on top of the original image. Make sure that the results are satisfactory. Print

out a suitable subset of the resulting images. Use this as a basis for your lab report and later exercises.

Note! As a refrence, your output images should look something like:

Question 7: Present your best results obtained with extractedge for house and tools.

## 6 Hough transform

In this exercise you will write a Matlab function houghline that uses the Hough transform to

determine linear approximations of a given number of curves, based on parameterizations of the type

x cos θ + y sin θ = ρ.

where x and y are suitably defined (e.g. centered) image coordinates, and ρ and θ have appropriately

selected definition areas with respect to the choices above. Your Matlab procedure should be declared

as

function [linepar acc] = …

houghline(curves, magnitude, …

nrho, ntheta, threshold, nlines, verbose)

where

• linepar is a list of (ρ, θ) parameters for each line segment,

• acc is the accumulator matrix of the Hough transform,

• curves are the polygons from which the transform is to be computed,

• magnitude is an image with one intensity value per pixel

(in exercise 6.2 you will here give the gradient magnitude as an argument),

• nrho is the number of accumulators in the ρ direction,

• nthetais the number of accumulators in the θ direction,

• threshold is the lowest value allowed for the given magnitude,

• nlines is the number of lines to be extracted,

• verbose denotes the degree of extra information and figures that will be shown.

Based on this function and extractedge you will then write a function houghedgeline that first

performs an edge detection step and then applies a Hough transform to the result. As a suggestion

this procedure might have a declaration similar to

function [linepar acc] = …

houghedgeline(pic, scale, gradmagnthreshold, …

nrho, ntheta, nlines, verbose)

where the (additional) arguments have the following meaning

• pic is the grey-level image,

• scale is the scale at which edges are detected,

• gradmagnthreshold is the threshold of the gradient magnitude.

To visualize the processing steps, you ought to show the important partial results on the screen.

Preferably, control the amount of partial results that is produced by varying the argument verbose of

the procedure, such that the value zero results in a “silent” procedure and the amount of information

displayed during execution increases with an increasing value of verbose. As an example of this, see

“type printcurves”.

To visualize the final result you should draw the lines overlaid on the original image. The easiest

way to do this is by creating polygons in the same format that was used by overlaycurves. If you

have enough time, add the possibility to draw only those line segments that correspond to points that

are close to the edge elements.

6.1 Hints and practical advice

• Think through the parametrization used for the edge segments (the choice of origin and definition

areas for ρ and θ).

Especially, make sure that you make use of the quantized accumulator space reasonably efficiently, and that discontinuities in the parametrization do not lead to any destructive effects.

• If you strike any programming problems, remember to divide the problem into smaller subproblems and test these partial modules separately. For example, you can be greatly helped by

writing test routines that draw a number of lines given values of ρ and θ.

• The interesting quantization of ρ corresponds to accumulator cells with ∆ρ in the neighborhood

of one pixel, while the interesting quantization of θ typically corresponds to accumulator cells

with ∆θ of about a degree.

• Large values of ∆θ considerably speed up the Hough transform. This fact is of value when

testing and debugging the procedure. Of the same reason, test and debug the function using

small images.

• You most easily detect local maxima in the accumulator by calling the function locmax8 in the

course library. The following code segment may be used to detect the local maxima and create

an index vector from these:

[pos value] = locmax8(acctmp);

[dummy indexvector] = sort(value);

nmaxima = size(value, 1);

Thereafter you can extract the index values in the accumulator that correspond to the nlines

strongest responses with a call like

for idx = 1:nlines

rhoidxacc = pos(indexvector(nmaxima – idx + 1), 1);

thetaidxacc = pos(indexvector(nmaxima – idx + 1), 2);

end

• Coordinate systems in curves versus images: When you want to go through the edge polygons,

the programming might be simplified if you start with the code in e.g. pixelplotcurves.m as

a template.

Regarding conventions for coordinate systems in images versus curves, you may look at the code

in pixelplotcurves.m (that accepts a given curve and writes it out as a collection of pixels).

function pixels = pixelplotcurves(image, curves, value)

insize = size(curves, 2);

trypointer = 1;

numcurves = 0;

while trypointer <= insize,

polylength = curves(2, trypointer);

numcurves = numcurves + 1;

trypointer = trypointer + 1;

for polyidx = 1:polylength

x = curves(2, trypointer);

y = curves(1, trypointer);

image(x, y) = value;

trypointer = trypointer + 1;

end

end

pixels = image;

Observe that the coordinate systems differ between images and curves. If we use the convention

that image elements are indexed as image(x,y), then the curve format has the x coordinate in

its second row and the y in the first.

The convention in the routines of the course library are such that the behaviors of following

operations are consistent:

t = triangle128;

showgrey(t)

c = zerocrosscurves(t-128);

overlaycurves(t, c)

showgrey(printcurves(t, c, -128));

Note that this definition of x and y axes is different from the notation used in “help contourc”.

• Quantization in accumulator space: For large values of ∆ρ and ∆θ the lines will be of low

accuracy. However, small values may result in multiple responses for the same line.

If you unable to achieve enough accuracy without getting too many multiple responses, it might

be advantageous to smoothen the histogram with a call to binsepsmoothiter.m, before detecting

local maxima.

• The output format of the results: The easiest way to visualize the results is by generating three

points of each extracted line. Note that Matlab automatically cuts off those parts of a line that

are located outside the image when the function overlaycurves is used, which means that you

do not explicitly have to compute where the line begins and ends. You can generate lines that

illustrate the results as follows:

x0 = < fill in >;

y0 = < fill in >;

dx = < fill in >;

dy = < fill in >;

outcurves(1, 4*(idx-1) + 1) = 0; % level, not significant

outcurves(2, 4*(idx-1) + 1) = 3; % number of points in the curve

outcurves(2, 4*(idx-1) + 2) = x0-dx;

outcurves(1, 4*(idx-1) + 2) = y0-dy;

outcurves(2, 4*(idx-1) + 3) = x0;

outcurves(1, 4*(idx-1) + 3) = y0;

outcurves(2, 4*(idx-1) + 4) = x0+dx;

outcurves(1, 4*(idx-1) + 4) = y0+dy;

In other words you have to convert each line defined by (ρ, θ) back to cartesian coordinates (x, y).

For convenience we decompose each line into 2 segments originated in (x0, y0), defined by 3 points

in a curve which forms a single line on the screen. You can adopt a different representation but

keep preferably the contourc format used by overlaycurves. Think of how you specified your

coordinate system in the Hough space!

• Test run the algorithm on small images (see help cutout). Exploit the interactivity of Matlab,

as well as the graphics facilities, if you need to debug your code.

Before you apply the algorithm to real images, make sure that it gives the correct results on

simpler test images, such as

testimage1 = triangle128;

smalltest1 = binsubsample(testimage1);

testimage2 = houghtest256;

smalltest2 = binsubsample(binsubsample(testimage2));

Apply your procedure houghedgeline to the simple test image triangle128.

Question 8: Identify the correspondences between the strongest peaks in the accumulator and

line segments in the output image. Doing so convince yourself that the implementation is correct.

Summarize the results in one or more figures.

Then apply the same procedure to the test image houghtest256. If the implementation is correct,

the results should be similar to those below:

Show the results of your procedure houghedgeline applied to the images few256, phonecalc256

and godthem256.

Original image Hough transform (10 lines)

Question 9: How do the results and computational time depend on the number of cells in the

accumulator?

6.2 Choice of accumulator incrementation function

A natural choice of accumulator increment is letting the increment always be equal to one. Alternatively, let the accumulator increment be proportional to a function of the gradient magnitude

∆S = h(|∇L|)

where |∇L| is the gradient magnitude and h is some monotonically increasing function. Furthermore,

you may in some cases use the local gradient direction.

Question 10: How do you propose to do this? Try out a function that you would suggest and see

if it improves the results. Does it?

6.3 A howto

When you develop function houghline, the follwing steps may help:

function [linepar, acc] = …

houghline(curves, magnitude, nrho, ntheta, threshold, nlines, verbose)

% Check if input appear to be valid

% Allocate accumulator space

% Define a coordinate system in the accumulator space

% Loop over all the input curves (cf. pixelplotcurves)

% For each point on each curve

% Check if valid point with respect to threshold

% Optionally, keep value from magnitude image

% Loop over a set of theta values

% Compute rho for each theta value

% Compute index values in the accumulator space

% Update the accumulator

% Extract local maxima from the accumulator

% Delimit the number of responses if necessary

% Compute a line for each one of the strongest responses in the accumulator

% Overlay these curves on the gradient magnitude image

% Return the output data