# Circle detection using invariance kernels

I'm not sure, if this belongs more to the signal processing forum or if I can post it in here, but I want to solve it via mathematica and it's some kind of image-processing I guess.

I have a noisy picture with a circle in it:

Now I want to detect the circle, by using so called invariance Kernels like a Fourier-Mellin transform. In short, I want to replicate what the authors of this paper (see Section 4) did in mathematica. So what they did, was using Fourier-Mellin transform to detect a noisy image with a circle in it. I have some troubles grasping the mathematics behind it and I hope someone could help me out in understanding and implementing it in a mathematica code.

• @holistic To get help for this kind of "where to even begin" type problem it should be original compared to already existing similar answers in this forum. If you little search this forum you will find many image feature detection related questions and plenty of them have excellent answers to help you with THE start. – PlatoManiac Feb 5 '14 at 13:10
• I think we need to be careful here. The reason why there was a problem with the link is because the paper is copyrighted by the publisher and not available for free. Does the OP have the rights to provide free access to this paper through M.SE and his/her dropbox account? – bobthechemist Feb 6 '14 at 0:06
• @bobthechemist I'm sure that the authors of the paper wouldn't mind, and I think there are many of us who are quite unhappy about the behaviour of scientific publishers. – Szabolcs Feb 6 '14 at 0:52
• @Szabolcs even the authors can't distribute on their own site. They transfer the publishing copyright to publisher. I think we are all unhappy except publishers. .. – s.s.o Feb 6 '14 at 1:10
• @s.s.o actually, reading through Elsevier's guidelines, I disagree. They can, with caveats, post on a website, just not the final form. – rcollyer Feb 6 '14 at 15:03

This is based on a reading of that paper. It's a ways from my areas and I do not claim to get it correct, but the code might be of use in any case.

Import the image.

im = Import["http://i.stack.imgur.com/POS2E.jpg"]


Get the basic data.

imat = ImageData[im];
idims = ImageDimensions[im];


Prepare a kernel matrix to convolve with. I start with a disk.

kmatinit = DiskMatrix[Min[idims]/2];
klen = Length[kmatinit];
posns = Flatten[Outer[List, Range[klen], Range[klen]], 1];
nf = Nearest[posns];


Now empty the center. Then I'll put concentric rings at different radii into different phases, that is, multiply by the jth ring by Exp[2*I*Pi*j/len] where len is the number of rings.

mid = Ceiling[klen/2];
minr = 5;
maxr = mid;
rrange = maxr - minr;
rvals = Exp[Range[0, rrange - 1]*2.*Pi*I/(rrange)];

kmat = kmatinit;
pts = nf[{mid, mid}, {Infinity, minr}];
Map[(kmat[[#[[1]], #[[2]]]] = 0) &, pts];
Do[
newpts = nf[{mid, mid}, {Infinity, minr + j}];
newpts = Complement[newpts, pts];
Map[(kmat[[#[[1]], #[[2]]]] *= rvals[[j]]) &, newpts];
pts = Join[pts, newpts];
, {j, rrange}];


Here is how it looks.

MatrixPlot[kmat]


Now convolve. If I understand correctly, the largest value should be in the vicinity of the circle's center (adjust due to convolution dimensions), and its phase should give us the radius.

lc = Chop[ListConvolve[kmat, imat, {1, -1}, 0]];
Dimensions[lc]
Max[Abs[lc]]
abslc = Max[Abs[lc]] - Abs[lc];
pos = Position[abslc, 0.]
p1 = pos[[1, 1]];
p2 = pos[[1, 2]];
lcpos = lc[[p1, p2]];
arg = lcpos/Abs[lcpos]

(* {200, 200}
949.436832261
{{97, 94}}
-0.336662959549 - 0.941625218262 I *)


So the center should be in the whereabouts of {97,94}/2. For the radius:

In[286]:= Nearest[rvals -> Range[Length[rvals]], arg]

Out[286]= {33}

It is the 33rd value in our set. We began at radius of 6 so this indicates our radius is 6+33-1 or 38. From the picture and the Manipulate from @Stephen Luttrell I believe it is actually around 40-41 or so, so this is not a terrible outcome. That said, I confess I don't really know what I'm doing.

--- edit ---

Here are some additional remarks, partly in response to a question.

(1) I make no claim to actually understanding that paper.

(2) The multiplications are so that concentric rings in the kernel array will have different phases. The idea is that rings (circles) that are not of the right radius and/or not centered with the circle we're seeking will tend to give noise that partly cancels. The ring of correct radius, when correctly centered, will all be in phase and not cancel with itself. We later use the phase information to figure out which radius ring this was.

It falls under the "phase coded annulus" description in the actual paper.

(3) There are probably flaws insofar as the different rings will be weighted according to their radii, so there may be a bias toward larger circles. I was not able to adjust for that in any way that worked well.

(4) If in doubt as to whether what I did makes sense, refer to remark (1). But rereading section 1.2.3 gives me the impression that at least I got the phase part correct.

--- end edit ---

• In your third code block, should kmatinit be kmat? kmat is not initialized. – Szabolcs Feb 7 '14 at 19:32
• @Szabolcs Yes. Now fixed. – Daniel Lichtblau Feb 7 '14 at 19:46
• Thanks @DanielLichtblau. Why did you multiply it with Exp[2*IPij/len]? To which part of the paper did you refer? – holistic Feb 10 '14 at 18:56
• I added some remarks that may help to explain.. – Daniel Lichtblau Feb 10 '14 at 19:13
• @DanielLichtblau you're not happy until you've sneaked in some Do, are ye? – Yves Klett Feb 10 '14 at 20:19