We'll use some operators from Chapter 3. If you want to write images out or to see them as pictures, you'll need to normalise them:

and we'll often need to clear them

From Chapter 4, we'll start using edge magnitude now, so we need the Sobel operator, and edge magnitude and direction

Chapter 5 Feature Extraction: CHAPTER5.MCD

We'll need some of the (Fourier transform) operators from Chapter 2:

Feature Extraction & Image Processing in Computer Vision

M. S. Nixon and A. S. Aguado

Academic Press, 3rd Edition 2012

This worksheet is the companion to Chapter 5 and implements the basic feature extraction techniques described therein. The worksheet follows the text directly.

First, we require a routine to match a template to an image, template matching. This is implemented by template convolution using a match functional. For each point within the template, if the point value matches the template value, then the votes are increased by one. The total of votes is then stored in an accumulator array which collects evidence. The point in the accumulator array where the votes are stored is the current position of the template. Accordingly, the accumulator array is two dimensional in terms of spatial position. This is then used to find the position of a shape in an image: the template is the shape and the maximum value in the accumulator array gives the position where the template best matched the image. Our operator, t_match, is then:

and sometimes we'll want to threshold them

initialise accumulator to zero

initialise range

address image

We'll need some noise

address template

increment votes if match

return accumulator

We'll need to check points are within an image

So let's have an image:

and a template of the rectangle we want to find in our image

let's have a look at them:

We'll need to find the co-ordinates of a point with a known value (brightness)

here's the image

and here's the template

So let's match the template:

and find the peak

and we'll need to shift an image around by

So we found our template at:

Let's look at the result, by putting a grey point on the location of the found shape:

and the image is:

The point of best match, the maximum in the accumulator space, is the starting point (upper left hand corner) of the template. That's why the position seems to be slightly out. It could have been plotted at the template's centre by adding half the size of the template to the co-ordinates of the maximum point.

let's have a look at the accumulator space:

Here, the peak is the point of the best match. The accumulator space shows that there was a some match elsewhere, between the template and the image. Note the large border space where no matching process occurs since the template would extend beyond the image.

We can also implement template matching by the Fourier transform. First, we need the template to be the same size as the image to which it is to be matched. We then need to pad out the template to be of the same size, inserting white at pixels where the template has no value (remember that our rectangle is black):

Let's pad....

Multiplication in the Fourier transform domain is actually equivalent to convolving two images. However, we want to correlate them to work out the point of best match. If we flip one of the images (reversing both its axes), then multiplication of the transform of the flipped image by the transform of the original image then gives correlation, template matching. So we need a routine to flip the template:

So let's have a peep at the flipped and padded template:

We'll be plotting Fourier transforms as their logarithm (to see what's in them) using

Original template

Padded template

Flipped and padded template

Now we can multiply the transform of the image by the transform of the flipped and padded template to implement template matching in the frequency domain by taking transforms

now we'll multiply them together:

and now we'll inverse Fourier transform to get the result

Note that the border regions are now occupied. A basic assumption of the DFT is that images replicate to infinity, so the votes in the border regions could be interpreted as incorrect.

Note that we've only used the rearrange function for display purposes. Let's join all this together (but without the rearrange function) as the routine t_corr

So let's apply it:

Did you notice any improvement in speed (over t_match)? You should have done. The template is 8x8, which is sufficient size for the FFT algorithm to give advantage. Also, Mathcad's FFT routines are compiled, and so frequency domain implementations in Mathcad can appear supremely fast, but this is artifice.

Let's find the maximum:

So we found our template at:

it's in the same place

Now, the size of the picture and the template are

Out of interest, let's read in bigger versions of the image and the template (the ones in the book)

By image domain calculation, this would require a massive computation. Let's do it by Fourier instead:

we can then get the result in Fig 5.4(g)

For the hardy amongst you, we'll do these bigger images by template matching in the image domain. It takes so long, that I switched it off. To make it go, select the operation (click on the rectangle after acc_space) and then press colon and off it'll go (as will you for a cuppa.....)

So far, we've only found the position of shapes. If we were required to find the shape irrespective of its size, then we would need to mach different sized templates to the image. In turn, this implies a three dimensional accumulator space, where the three dimensions are the spatial co-ordinates and the size. If we were to progress to orientation, to find a shape irrespective of its rotation, then we would increase the dimensions of the accumulator space by one more. Perspective transformation would increase it even more. These are only modifications to the basic scheme: they do not change its basis, only its implementation.

We'll now see what happens with noisy pictures, let's read in a noisy version of our original image

now let's match the same template to it:

and compare the match result with that for the original (clean) image, first the peaks:

so the peak is reduced (as you expect). Now let's look at the location:

but it's still in the right place. This is one of the virtues of template matching. It finds shapes in the right place, even when there's a lot of noise. Let's see it again, for some varaitions in the noise level

Drawing a black cross on a white background where the origin of the shape is found, we have

Good, all fine and dandy

For the more noisy image

Still working, but hard for human vision here.

So it's still OK. Let's try the noisiest image:

OK, we've broken it- which is what we'd expect for such a noisy image. In fact, we might not have -it still wolrks occasionally. Just click on the equation generating the noisy image (noisy_picture6:=) and press F9. This gives a diferent instantiation of the noise. Occasionally it works, but not very often.

Further, we can still detect the shape when it is occluded (hidden) by another shape.For this, we shall use a function which draws a vertical grey bar on the image as

So let's draw the bar, in the first case, not occluding the shape

So we'll increase the level of occlusion

It's partly hidden now

Ah, still working

So we'll increase the level of occlusion even further

It's fully hidden now

Ah, broken now - as expected

The good performance is because it is a version of matched filtering, delivering a result which is optimal by the least squares criterion. We could make the implementation more robust by using edge detection (thereby concentrating on a shape's borders), but we can only make it faster by Fourier implementation.

A template isn't the only shape parameterisation. The Histogram of Oriented Gradients operator describes shapes by using local properties of edge direction, in fact rather as its name implies. It uses basic edge detection

For an image

We derive edges

and we measure edge directions as

Giving, in image form

If we measure these directions in 8 bins

So we now have a descriptor we can match to the image

The other way to make template matching faster is to notice that it is actually a mapping transform. This is the Hough Transform (HT). We'll start the HT analysis by looking at the HT for lines.

The basis of the HT for lines is the duality between a line in x,y co-ordinate space, and a line in m,c co-ordinate space. The equation for a line is:

y = mx + c

for points with co-ordinates x,y collinear with gradient m and intercepting the y axis for a value c. If we change this round we get

c = -xm + y

for points with co-ordinates c,m collinear with gradient -x and intercepting the c axis for a value y. So a point x,y maps to a line in m,c. The co-ordinate space m,c is the accumulator space; lines drawn in the accumulator space are lines of votes. So the co-ordinates of a point x,y give, for all possible values of m, a set of possible values for c. The location of each accumulator array point for all points m,c is incremented for each vote.

We'll need some lines to test our algorithms on. Let's have a subroutine which generates lines

and we'll accumulate evidence via the HT for lines A_line then processes the edge image E_mag according to :

clear accumulator space

threshold for edge magnitude

look at whole image

if edge point is above threshold

for all possible values of m

calculate value of c

and increment accumulator cell

dummy operation for clear layout

return accumulator

Let's synthesise an edge image as a perfect line of brightness value 255 by calling the line function

This is the image of our line. It contains 36 points.

and let's generate the votes for it

and see them

Note that the peak is the same as the number of points in the original line

The maximum in the accumulator space is

These are the values which we used to draw the original line. So it works! Try another line (by using different values for m and c), and check that it still works.

The major difficulty with the basic HT for lines is in definition of an appropriate parameter space. The one just used uses vales of m between 0 and the image size to calculate values of c. In fact, m can be infinite for a vertical line. This would imply a rather large parameter space. To avoid it we use the foot-of-normal parameterisation where lines are represented by

xcosq + ysinq = r

where x,y are the co-ordinates in image space and r,q are the co-ordinates of the accumulator space. This equation is used to replace the voting functional c = -xm + y. This gives A_polar

set resolution for q

find maximum r

clear accumulator

set edge threshold

cover whole image

if edge exceeds threshold

for all values of q

calculate equivalent radians

calculate r

increment accumulator cell if co-ordinates o.k.

dummy operation for legible layout

return accumulator

So let's try it out:

So the relevant bit of the accumulator now looks like

We'd better check this peak out! How can we tell if it's right (especially when there's another pretty close peak). We need a routine to draw lines by using polar format as

Original line

Result from the HT

The result is fine. It's just that our line straddles two adjacent cells in the accumulator space. The solution is to make the quantisation of the accumulator space smaller

So let's move on to the HT for circles. This uses the same duality as the HT for lines: a circle of points with co-ordinates x,y centre x 0,y 0 radius r is

(x - x 0)2 + (y - y 0)2 = r2

which is also a circle with co-ordinates x 0,y 0 centre x,y radius r. So the equation again gives a voting functional, but this time in a three dimensional accumulator space, where the three dimensions are position (x,y) and the radius. Unfortunately, 3D matrices are not routine in Mathcad. We'll first need a routine to generate circles as:

We'll need a 64x64 image, let's read in the eye image to declare one

So, for a fixed radius we require a 2D accumulator space (equivalent to a slice through the 3D space) giving the routine A_circle as

clear accumulator

cover picture

if edge point exceeds threshold

use angular increments

to generate a circle of votes,

the accumulator cells to be incremented

dummy operation for legibility

and return the accumulator

So let's try it at the radius we specified for the original circle

let's find the peak:

perfect! We'll chop off some of the accumulator array to save time when we see it:

Now let's look at the structure of the votes around the peak:

Now try different values of the radius - what happens to the peak? Also, have a look at the structure of the votes in the background, is it what you would expect?

Let's try it on a real image, another image of an eye:

now edge detect it

and threshold it

eye

now vote, trying a radius of 12

find the peak

thresholded edges

draw a circle at the peak

Try other values for the threshold and the radius. Can you explain your results? How could you improve on the results you achieve?

located circle

We could now move on to the HT for ellipses. These are defined as

(x - x 0)2 / a2 + (y - y 0)2 / b2 = 1

which is an ellipse with co-ordinates x,y centre x 0,y 0 with major and minor axes a and b, respectively. This implies that we need a four dimensional accumulator space; if we include rotation we need a five dimensional accumulator space. That could be rather large. It's only a minor modification to the HT for circles anyway - the voting functional just draws ellipses of votes, rather than circles. So we'll now move on to a technique which is not restricted to conic sections such as lines, circles or ellipses.

The generalised HT can be used to find arbitrary shapes, shapes which are not expressed in closed analytic form, like lines and circles. Essentially it stores the voting functional as a look up table. The look up table maps image points to the accumulator space and finding a shape is then equivalent to finding the peak in the accumulator array. The look-up table depends on edge direction. The edge direction at a point of the template shape indexes the distance to, and direction of, a reference point. When analysing a target shape, the edge direction at a point is used to vote for possible centre points. Consider finding a shape like a card-suit spade in an image. The template of the spade is:

We'll edge detect it using Sobel again

and threshold it

We now need to make the look-up table. This is actually called the R-table. It needs a reference point, we'll use a function to calculate the centre co-ordinates of our template shape.

initialise centre co-ordinates

address all points in template

if a point exists

tot up x total

and y total

and no. of points

return average as centre

We'll initialise the accumulator array cells to zero, the number of bins is the number of possible values of edge direction

Now we'll make the R-table

start with a vector of zeroed bins

get the reference point

for legibility later

address whole image of template

if template point

get distance to reference point

get direction to reference point

integerise edge direction

work out where to put new point

include it as direction and distance

insert new (next) zeroed element

and we invoke this as

and we accumulate evidence as

clear the accumulator

address the image

if edge point above threshold

evaluate integerised direction

for appropriate R-table point

evaluate distance

and direction

of point to be voted for

and vote

more cells?

Let's apply it. We'll add some noise to our template, and then shift and then try and find it again. First we'll add some noise to a shifted image:

Now we'll edge detect it

and threshold the shifted version

So our target image is now:

Could be tricky....let's see what the GHT can do.

and we'll find the peak

impressive eh! But remember, it's a result (nearly) equivalent to template matching. To run it again, with a different selection of noise points, just select the function noisy_template := addcondiments() and run it again. Each time, you'lll get a different pattern of noise. Check the extraction still works.

We better check with the original template:

and if we add our shift, we get the right value

This is pretty close. Out of interest, let's see the result by template matching. We're matching 31x31 images, so we'll use Fourier. We'll need to pad out our image (for the FFT) to 32x32 so we'll need pointers

and clear an image

now we'll copy our target image into it

This is the co-ordinates of the top left hand corner. The centre is at 10,10 so this again gives the right answer. The difference between template matching and the GHT is that the GHT uses edge direction, and a discrete set of bins. If the Fourier implementation of template matching was not artificially fast, you'd see no speed advantage by template matching.

Finally, we'll look at the votes around the peak

Not such a bad result, the adjacent peak is not quite so high, compared with the peak. Look at a printout to see the dispersion.

So we can now find arbitrary shapes too, in a manner providing a result which is equivalent to template matching, but faster. This is where we will leave evidence gathering techniques for the moment, and return to them later. We'll now move to techniques which evolve to the solution, rather than find it by post-processing an accumulator array.