Accuracy on detected radii, minor axes and major axes

[tobbelobb]

I'm experimenting with using ED_Lib for a motion tracking system that tracks colored spheres on images.

Today I conducted the following experiment:

1. Create a render where 135 red spheres of known size are placed out in a grid pattern on a plane, against a white background.
2. Use EDColor to detect as many as possible (never mind false positives, they are not a problem).
3. Use center point and minor axis (not major axis) gotten from EDColor to calculate a 3d position for each sphere.
4. For each position, place out a different colored, equal sized sphere alongside the red one.
5. Create a new render, which shows the error in 3d position, for each sphere.

The output of EDColor (result of step 2) was like this:

The image above is scaled down. The real output from step 2 is very large (123Mpix).

We see that EDColor regards spheres that are projected close to the image center, as well as halved sphere projections along the top and bottom edges of the image as circles. We also see that complete sphere projections far from the image center is recognized as ellipses.

So ED_Lib deviates a bit from reality. All the projections of spheres, except the middle one, are really ellipse shapes. This deviation was expected.

The result of step 5 was like this: or, seen from straight above:

Here we see that ED_Lib is good enough at finding the correct ellipse center. We also see that ED_Lib is good enough at finding the radii/minor axes of the projections that are close to the image center (nearly perfectly circular ones).

However, it looks like ED_Lib has troubles detecting accurate enough minor axis size for the more elliptical projections, both those that are approximated to be circles, and those who are recognized as ellipses.

In reality, the radii/minor axes of all the projections are identical.

Here is a close up of one of the most elliptical sphere projections, overlaid with the ED/EDPF edge image in black, and the detected ellipse, in yellow: We see that the detected ellipse is too small by ~1 pixel (some places 0, some places 2, but never too large). This causes the calculation in step 3 to find a position that is too far away from the camera.

Here is another example, with one of the most elliptical projections that is still recognized as a circle by EDCircle: We see a white gap between the green line and the red sphere, which means that the radius found by EDCircle is bigger than the real minor axis. This causes the calculation in step 3 to find a position that is too close to the camera.

I have thoroughly tested my equations, and done the same experiment with another circle detector to make sure the 3d position errors you see above really originate from ED_Lib's axis detection errors.

So finally, the question: I want to minimize or compensate for the radius/minor axis error. I can imagine a few ad-hoc solutions that would scale the radii or minor/major axes after the fact, but I wonder: Is there anything we can do to improve the accuracy of ED_Lib itself? Where would you start?

Best regards, and thanks again for ED_Lib! tobben

PS! I also did the experiment with step 3. changed to include the detected major axis in the calculations. That produced a slightly worse result, errors in the same direction, but larger, so never mind the major axis, I just need accurate enough minor axes.

PPS! My code is here. gitlab.com/tobben/hpm and the usage of it is here gitlab.com/tobben/hp-mark

[CihanTopal]

First of all there is something that confused me. Do you mean EDCircles by EDColor, because EDColor only returns edge segments, not higher level of information such as center or semi-major/minor axes? Or you detect edges and use them as the primary features to build up your own ellipses?

I found this problem really weird since ED extracts edge segments by literally picking edgels one by one with respect to gradient response of image details. Here is a demo video how it works https://youtu.be/-Bpb_OLfOts.

In another example, in this video (https://youtu.be/cPLmmVLGrdQ), we are building up pupil contour's based on ED edge segments and it works very well.

Similarly, we utilize ED outcome in STag marker and we showed that the localization of the edges is superior by comprehensive experimentation that you can see in the paper https://www.sciencedirect.com/science/article/pii/S0262885619300903.

So, currently I have no idea why this is happening. There might be some implementation discrepancies while we were upgrading the first C implementation to C++ for this github repo, but I don't think this is likely very much.

Is this bad localization problem happening only for ellipses, or it is the same for circles as well? If it is, then perhaps the issue is related to the ellipse fitting method that we utilized in ellipse fitting? Or, may be, EDColor performs the same edge segment extraction algorithm on the Lab color space, therefore, intensity overlap by projecting of RGB values onto the gray scale (3D to 1D) becomes less problematic. So, can those mislocalized edges are actually happen due to the inefficient gradient magnitude computation of your images in Lab color space? (I am just thinking aloud, I do not think this is likely either).

[tobbelobb]

Cool, I'm happy to hear that it's supposed to work differently. Yes, I'm also puzzled by the results, but it's a bit exciting and interesting as well.

Yes, I create an EDColor object from which I create an EDCircles object.

The edgeImage within the EDColor object (both the one built by the ED constructor, and the (filtered) one built by the EDColor constructor) places the ellipse perimeters exactly right.

The segments within the EDColor object also follows the edges of the edgeImage perfectly.

So the discrepancy emerges somewhere inside EDCircles.cpp, in the constructor that takes an EDColor as an argument.

Is this bad localization problem happening only for ellipses, or it is the same for circles as well?

The closer to circular the projection gets, the smaller the error gets, although it does not go completely to zero. Here is an example of the detection of a perfect circle:

We see that EDCircle has chosen an ever so slightly too small radius, even for the perfect circle. The black perimeter edge there is the edgeImage/segment input that EDCircle has used to calculate the yellow circle output.

So there is some logic inside EDCircle.cpp logic, both the circle fitting logic and the ellipse fitting logic, that is resposible for the ~1px discrepancy.

I understand that the algorithm calculates radii and centers of circle segments independently, and then tries to join circle segments that have similar radii and centers.

I wonder if afterwards, when a circle is determined to be found, is the final radius refined, based on all the included segments?

I have looked into EDCircle::CircleFit in some detail. I noticed that the (x,y) pairs of input, that we fit the circle to, contain integers, so there is a truncation error there. I did a hasty experiment where I moved the center, and also the truncation error of each pixel outwards from the circle radius, like so:

  for (int i = 0; i < N; i++) {
    xAvg += x[i] - 0.25;
    yAvg += y[i] - 0.25;
  }

 ...
  for (int i = 0; i < N; i++) {
    double u = ((x[i] < xAvg) ? x[i] - 0.5 : x[i] + 0.5) - xAvg;
    double v = ((y[i] < yAvg) ? y[i] - 0.5 : y[i] + 0.5) - yAvg;

... and I managed to get a much better, yet still not error free circle detection:

We see that the radius is now correct, but the center is still 1 or 0.5 px too close to the bottom right corner of the image.

[tobbelobb]

I should also mention that since the image is huge (123Mpx), the truncation errors of the floating point arithmetic within EDCircles::CircleFit get larger. They might play a role here. We should do an experiment with long double arithmetic.

But I'm off for today, talk to you tomorrow ;)

[tobbelobb]

I did another experiment today, where I compensated for the truncation error in EDCircle::CircleFit that I mentioned earlier. The center of each pixel (x, y) is actually at (x+0.5 px, y+0.5 px). Accounting for this lowered the center point error from 0.062 mm to 0.007 mm, ie almost by a full order of magnitude :)

In another experiment, I found that the floating point arithmetic truncation error is insignificant in all my test cases.

When it comes to the radius of the perfect circle, that can be fine tuned by adding or subtracting to the line

*pr = R + something;

For some reason that I don't fully understand, I get the best results when setting something = 1/3.

After these adjustments, the circle fit (in yellow) matches the edgeImage (in black) almost perfectly: Nowhere do we see red pixels in between the yellow and the black pixels anymore :tada:

With these adjustments, perfect circles are determined perfectly enough for my application. However, I will always see ellipses in my real use case. So I tried returning a high error from CircleFit regardless of the input. That way, the EllipseFit function will take over and fit ellipses everywhere. This gave a very positive result:

We see that the error in the minor axis is a constant. So I will dive into EDCircles::EllipseFit, and try to fine tune it.

I will have to turn CircleFit back on in the future though, since EllipseFit gives a lot more false negatives on real world images.

So, EDCircles is like two detectors in a sequence. The first detector finds a lot of circles, and is very rough for slightly elliptical circles, but with very few false negatives. The second detector finds ellipses, with a very predictable error, so it's not rough at all, it's very precise, but it gives a lot of false negatives.

An idea to combine the two would be if we could change:

if (isCircle) CircleFit();
else if (isEllipse) EllipseFit();

into

if (isCircle or isEllipse) EllipseFit();

This comment is too long now, sorry for that.

Cheers

[tobbelobb]

EllipseFit had the same trucation error. Got the same 0.062 -> 0.007 mm error reduction by adding the + 0.5 to these two lines in EllipseFit:

    tx = x[i - 1] + 0.5;
    ty = y[i - 1] + 0.5;

[tobbelobb]

And after also adding the + 1.0 to these lines inside EDCircles::ComputeEllipseCenterAndAxisLengths()

  // semimajor axis
  a = sqrt(F3 / A2) + 1.0;
  // semiminor axis
  b = sqrt(F3 / C2) + 1.0;

... we get this magnificent result:

We see that all the spheres that were not cut by the edge of the original image have been positioned back correctly. The error on the Z-axis is reduced from 12.99 mm to 0.52 mm. That is a 23x improvement for my usecase :tada:

[tobbelobb]

The mystical 0.33 value mentioned before has roots both in how the edge drawing in ED works, and how the renderer that created my test image works, as well as a definition question about how large part of the outermost colored pixels really belong to the circle.

I made a more controlled test, that removes the renderer from the equation, here: https://gitlab.com/tobben/hpm/-/blob/05d6f532fa1aa5f5afb643a7dea0bef0926d55ea/hpm/ed/ED.test.c++#L98

It turns out that whether to put R+0, R+0.33, or R+0.5 is hard to determine when I'm training on rasterized images. Because of the rasterization, some red will be outside of the circle, or some white will be inside the circle, or both. I've opted for R+0.5 in my code and made sure that test.cpp finds the exact same number of segments, circles, etc, as before.

The test is not perfect, but it is ok, and it would have caught the previous truncation error inside CircleFit.

Cheers

[CihanTopal]

I am sorry I couldn't catch up with your progress, but I will try to provide a summed up reply for your comments.

First,

if (isCircle or isEllipse) EllipseFit();

is not a good idea since ellipse fitting is way more computationally expensive than circle fitting.

Regarding integers in the following part,

for (int i = 0; i < N; i++) {
    xAvg += x[i] - 0.25;
    yAvg += y[i] - 0.25;
}

I checked it and they are all double. Are we talking about the same lines of code (EDCircles.cpp - 2937:2940) ?

May all these inaccuracies happen due to the fact that the pixels locations are actually meant by their top-left corner, where as you shifted that point to the physical center of the pixel by adding [0.5, 0.5] to them?

I am glad that you eventually made it worked for your problem!

[tobbelobb]

Hi!

May all these inaccuracies happen due to the fact that the pixels locations are actually meant by their top-left corner, where as you shifted that point to the physical center of the pixel by adding [0.5, 0.5] to them?

Exactly :) It doesn't solve all the inaccuracies, but some of them. The xy values are doubles, but they take on values that are very close to whole numbers, like 1.000, 560.000, and so on. Adding [0.5, 0.5], solved a big part of the issue.

is not a good idea since ellipse fitting is way more computationally expensive than circle fitting.

Yes, I agree that the EllipseFit() function is quite computationally expensive. I could use a cheaper, more ad-hoc compensation approach on the application side (not inside ED_Lib), I guess.

I think it would also be feasible to refine a circle fit into an ellipse fit quite cheaply. We could fixate the center and, since the circle was found by least squares fitting, we know that the found radius is roughly equally far away from the real minor and the major axes. So we could do a couple of Newton-Raphson iterations, searching for only one variable x.

semiMinor = radius - x;
semiMajor = radius + x;

We know that x is very small, and we could use the fitting error from CircleFit to create a very good first guess. Maybe we could just use the guess, and don't need any Newton-Raphson iterations at all. I'll try it.

Thanks for your feedback!

[tobbelobb]

Was able to find x.

So, I store the error from CircleFit as a member in the mCircle class. Then I convert a circle into an ellipse like this in my application code:

  Ellipse(mCircle const &circle)
      : m_center(circle.center),
        ...,
        m_minor(2.0 * (circle.r - (sqrt(3) * circle.err) + 0.5)), m_rot(0.0) {}

The sqrt(3) converts the root mean square calculated inside CircleFit into the mean radius error. Integrate x² from 0 to x' and divide by x', and take the square root to get that sqrt(3) analytically. And as always, a ~0.5 term representing truncation sneaks in.

Using similar logic inside ED_Lib directly might or might not make sense, I don't know.

My final victory picture:

I hope it has been fun for you as well. I will close this issue now. Reopen if you have comments, thoughts, or questions.

BR