Smooth a bumpy circle

Question

I am detecting edges of round objects and am obtaining "bumpy" irregular edges. Is there away to smooth the edges so that I have a more uniform shape?

For example, in the code below I generate a "bumpy" circle (left). Is there a smoothing or moving average kind of function I could use to obtain or approximate the "smooth" circle (right). Preferably with some sort of parameter I can control as my actual images arn't perfectly circular.

import numpy as np
import matplotlib.pyplot as plt

fig, (bumpy, smooth) = plt.subplots(ncols=2, figsize=(14, 7))

an = np.linspace(0, 2 * np.pi, 100) 

bumpy.plot(3 * np.cos(an) + np.random.normal(0,.03,100), 3 * np.sin(an) + np.random.normal(0,.03,100))

smooth.plot(3 * np.cos(an), 3 * np.sin(an))

Answer 1

If the shapes in question approximate ellipses, and you’d like to force them to be actual ellipses, you can easily fit an ellipse by computing the moments of inertia of the set of points, and figure the parameters of the ellipse from those. This Q&A shows how to accomplish that using MATLAB, the code is easy to translate to Python:

import numpy as np
import matplotlib.pyplot as plt

# Some test data (from Question, but with offset added):
an = np.linspace(0, 2 * np.pi, 100) 
x = 3 * np.cos(an) + np.random.normal(0,.03,100) + 3.8
y = 3 * np.sin(an) + np.random.normal(0,.03,100) + 5.4

plt.plot(x, y)

# Approximate the ellipse:
L, V = np.linalg.eig(np.cov(x, y))
r = np.sqrt(2*L)                       # radius
cos_phi = V[0, 0]
sin_phi = V[1, 0]                      # two components of minor axis direction
m = np.array([np.mean(x), np.mean(y)]) # centroid

# Draw the ellipse:
x_approx = r[0] * np.cos(an) * cos_phi - r[1] * np.sin(an) * sin_phi + m[0];
y_approx = r[0] * np.cos(an) * sin_phi + r[1] * np.sin(an) * cos_phi + m[1];

plt.plot(x_approx, y_approx, 'r')
plt.show()

The centroid calculation above works because the points are uniformly distributed around the ellipse. If this is not the case, a slightly more complex centroid computation is needed.

Answer 2

You can do this in frequency domain. Take the (x,y) coordinates of the points of your curve and construct the signal as signal = x + yj, then take the Fourier transform of this signal. Filter out the high frequency components, then take the inverse Fourier transform and you'll get a smooth curve. You can control the smoothness by adjusting the cutoff frequency.

Here's an example:

import numpy as np
from matplotlib import pyplot as plt

r = 3
theta = np.linspace(0, 2 * np.pi, 100) 
noise_level = 2
# construct the signal
x = r *  np.cos(theta) + noise_level * np.random.normal(0,.03,100)
y = r *  np.sin(theta) + noise_level * np.random.normal(0,.03,100)
signal = x + 1j*y
# FFT and frequencies
fft = np.fft.fft(signal)
freq = np.fft.fftfreq(signal.shape[-1])
# filter
cutoff = 0.1
fft[np.abs(freq) > cutoff] = 0
# IFFT
signal_filt = np.fft.ifft(fft)

plt.figure()
plt.subplot(121)
plt.axis('equal')
plt.plot(x, y, label='Noisy')
plt.subplot(122)
plt.axis('equal')
plt.plot(signal_filt.real, signal_filt.imag, label='Smooth')

Answer 3

I recommend to use FIR filter. The idea behind is to use weighted average of points arround filtered point... You just do something like

p(i) = 0.5*p(i) + 0.25*p(i-1) + 0.25*p(i+1) 

where p(i) is i-th point. Its a good idea to remember original p(i) as it is used for next iteration (to avoid shifting). You can do each axis separately. You can use any weights but their sum should be 1.0. You can use any number of neighbors not just 2 (but in such case you need to remember more points). Symmetric weights will reduce shifting. You can apply FIR multiple times ...

Here small 2D C++ example:

//---------------------------------------------------------------------------
const int n=50;     // points
float pnt[n][2];    // points x,y ...
//---------------------------------------------------------------------------
void pnt_init()
    {
    int i;
    float a,da=2.0*M_PI/float(n),r;
    Randomize();
    for (a=0.0,i=0;i<n;i++,a+=da)
        {
        r=0.75+(0.2*Random());
        pnt[i][0]=r*cos(a);
        pnt[i][1]=r*sin(a);
        }
    }
//---------------------------------------------------------------------------
void pnt_smooth()
    {
    int i,j;
    float p0[2],*p1,*p2,tmp[2],aa0[2],aa1[2],bb0[2],bb1[2];
    // bb = original BBOX
    for (j=0;j<2;j++) { bb0[j]=pnt[0][j]; bb1[j]=pnt[0][j]; }
    for (i=0;i<n;i++)
     for (p1=pnt[i],j=0;j<2;j++)
        {
        if (bb0[j]>p1[j]) bb0[j]=p1[j];
        if (bb1[j]<p1[j]) bb1[j]=p1[j];
        }
    // FIR filter
    for (j=0;j<2;j++) p0[j]=pnt[n-1][j];                    // remember p[i-1]
    p1=pnt[0]; p2=pnt[1];                                   // pointers to p[i],p[i+1]
    for (i=0;i<n;i++,p1=p2,p2=pnt[(i+1)%n])
        {
        for (j=0;j<2;j++)
            {
            tmp[j]=p1[j];                                   // store original p[i]
            p1[j]=(0.1*p0[j]) + (0.8*p1[j]) + (0.1*p2[j]);  // p[i] = FIR(p0,p1,p2)
            p0[j]=tmp[j];                                   // remeber original p1 as p[i-1] for next iteration
            }
        }
    // aa = new BBOX
    for (j=0;j<2;j++) { aa0[j]=pnt[0][j]; aa1[j]=pnt[0][j]; }
    for (i=0;i<n;i++)
     for (p1=pnt[i],j=0;j<2;j++)
        {
        if (aa0[j]>p1[j]) aa0[j]=p1[j];
        if (aa1[j]<p1[j]) aa1[j]=p1[j];
        }
    // compute scale transform aa -> bb
    for (j=0;j<2;j++) tmp[j]=(bb1[j]-bb0[j])/(aa1[j]-aa0[j]);   // scale
    // convert aa -> bb
    for (i=0;i<n;i++)
     for (p1=pnt[i],j=0;j<2;j++)
      p1[j]=bb0[0]+((p1[j]-aa0[j])*tmp[j]);
    }
//---------------------------------------------------------------------------

I added also checking BBOX before and after smoothing so the shape does not change size and position. In some cases centroid is better than BBOX for the position correction.

Here preview of multiple application of FIR filter: