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