Circle detection using invariance kernels

Question

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.

Answer 1

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 ---