How to apply a fit on this shape?

Since my data doens't really count as a function cftool is not an option. I wonder whether this is still faisable to fit:

(I tried several elliptical fits but they don't seem to work. May be because I dealt with it a little clumsily)

Answer 1

[John D'Errico]

What does it mean to "fit" it to you? What will you do with that fit? The data is not a function, just a set of points, that happen to represent a closed curve in the plane.

A simple solution might be to convert to polar coordinates.

help cart2pol

Now, fit the function r(theta). Be careful, since r is a periodic function. It should have the property that r(0) == r(2*pi). So a polynomial fit would be a bad idea. Instead, perhaps fit this as a sum.

r(theta) = a0 + a1*sin(theta) + b1*cos(theta) + a2*sin(2*theta) + b2*cos(2*theta) + ...

Essentially, I am describing what is essentially a Fourier series there (that is, not an fft, but a Fourier seies.)

https://en.wikipedia.org/wiki/Fourier_series

As it is constructed, the function ALWAYS has the property of being periodic. And, while you could, in theory, trivially use the curve fitting toolbox for this task, there is no reason to do so. All of the parameters can be estimated using one call to backslash.

x = x(:);
y = y(:);
mux = mean(x);
muy = mean(y);
[Theta,R] = cart2pol(x- mux,y-muy);
[Theta,ind] = sort(Theta);
R = R(ind);

plot(Theta,R,'o')

n = numel(x);
M = 6;
A = [ones(n,1),sin(Theta*(1:M)),cos(Theta*(1:M))];
coeffs = A\R;

Rhat = A*coeffs;

hold on
plot(Theta,Rhat)

Larger values of M will yield a better approximation. A perfect circle would have m == 0. Because your curve has an oval shape, you need to have a few extra terms, but not many. M == 6 will require the estimation of 13 coefficients for the series.

I lack your data, so I generated some random crap data as an example, using ginput.

plot(x,y,'o')

Nothing exciting. Vaguely ovoid. I'm actually surprised it was that good.

R as a function of Theta seems reasonable.

And now, we can reconstruct the original curve.

[xhat,yhat] = pol2cart(Theta,Rhat);
plot(x,y,'bo',xhat + mux,yhat + muy,'-xr')

Surprisingly good for a simple, moderately low order Fourier series approximation, given the data came from my shaky eye/mouse coordinated hand.

[Niklas Kurz]

Big up, that's some really nice representation. I'm truely astonished. Unfortunaly I lack of deep mathimatical knowledge (At the moment I'm just being in 2nd semester, Fourier Transformation will follow at least in the next one).

However, I attached my data above. I already tried your code receiving a really neat figure. In contrast to other elliptical fits in add-ons.

So to cut a long story short: Actually u answered my question excellently, but I just want to calculate the surface inside this shape. Therefore I thought I'd be smart to draw an ellipse function out of my Data and simply use the parameters for the known area of an elilpse.

But this might be inconvenient now. So how u'd solve it instead posessing such a great toolbox of mathematical skills?

[John D'Errico]

Are you asking to compute the area contained in the region found from the fit?

function [a0,sincoefs,coscoefs,mux,muy,Rfun,area] = FSfit(x,y,nterms,plotflag)
  % approximates the closed curve in (x,y) using a sine & cosine series
  % the number of sine & cosine terms is indicated by nterms.
  % if plotflag is supplied as true then plots are generated too
  
  if (nargin < 4) || isempty(plotflag)
    plotflag = false;
  end
  
  x = x(:);
  y = y(:);
  nx = numel(x);
  
  mux = mean(x);
  muy = mean(y);
  [Theta,R] = cart2pol(x- mux,y-muy);
  
  A = [ones(nx,1),sin(Theta*(1:nterms)),cos(Theta*(1:nterms))];
  
  coefs = A\R;
  
  a0 = coefs(1);
  sincoefs = coefs(2:nterms+1);
  coscoefs = coefs(nterms+2:end);
  
  % the area contained in the region is actually easy to compute
  area = pi*(a0^2 + sum(sincoefs.^2)/2 + sum(coscoefs.^2)/2);
  
  Rfun = @(th) reshape(a0 + sin(th(:)*(1:nterms))*sincoefs + cos(th(:)*(1:nterms))*coscoefs,size(th));
  
  if plotflag
    th = linspace(-pi,pi,500)';
    figure
    plot(Theta,R,'b.',th,Rfun(th),'r-')
    
    figure
    plot(x,y,'b.',Rfun(th).*cos(th) + mux,Rfun(th).*sin(th) + muy,'r-')
    axis equal
  end

After trying smaller values, it looks like nterms as 6 seems a pretty good approximation to the data.

[a0,sincoef,coscoef,mux,muy,Rfun,area] = FSfit(x,y,6,true);

The area of this region is easy enough to compute. I'll use the symbolic toolbox for no good reason except I'm too lazy to think. We want to form the integral r(theta)*r(theta)*dtheta, where theta varies from -pi to pi.

syms t a0 a1 a2 b1 b2
R = a0 + a1*sin(t) + a2*sin(2*t) + b1*cos(t) + b2*cos(2*t);
int(R^2/2,[-pi,pi])
ans =
pi*(2*a0^2 + a1^2 + a2^2 + b1^2 + b2^2)/2

Do you see what happens? We have a term in there that looks like the area of a circle, and then we add in terms with the squares of each trig coefficient. That actually makes mathematical sense. Remember that sin and cos are orthogonal functions on [-pi,pi]. In fact sin(theta) is also orthogonal to sin(k*theta) for integer k values of k>1). Almost everything drops out. Ok, I said remember, but if you've not done anything with Fourier series yet, the orthogonal function stuff may be new to you. Such is life.

And just to show the result is exact, we can compute it both using the formula, and using integral. Note that to compute the area, mux and muy never factor in.

area
area =
          2.92889047909024
          
integral(@(t) Rfun(t).^2/2,-pi,pi)
ans =
          2.92889047909024

That looks about right for the area too, a little less than pi.

[Niklas Kurz]

Wow, thank you so much for this explanation and your friendly advices, it really helps me deepening my knowledge. :) All the best for you too!

[Niklas Kurz]

If I may ask additionally: How to determine the area beneath the curve? Pretending the ellipse would be above the x axis like in the data I attached as well. Since the values are repeating trapz(X,Y) doesn't seem right.

[John D'Errico]

Why does this question seem to be veering in random directions, with me trying to chase a moving target?

If the goal is to find the area betwen the lower lobe of that closed curve, and the x-axis, there is a problem. You need to find the right and left hand end points to integrate between. A problem is the data is relatively noisy, but also, the noise is not really true noise.

And given this data, the code I wrote needs quite a few terms to perform even reasonably well.

FSfit(x,y,100,1);

While the fit does not look too bad when viewed as an entire plot, we can see it shows serious oscillations in places. Another option is to use a smoothing spline, but that looks just as bad.

I'd suggest a big part of the problem is the data is heavily quantized in the x variable. For example, there are 9287 data points in your curve. However, x is clearly not any sort of continuous variable.

numel(unique(x))
ans =
   125
   
numel(unique(y))
ans =
        9287

One thing that I often do whenever someone asks a question about how to model some data, the first thing I often ask is what do you want to get out of it? I should have pushed for that here. If the goal is merely to find the area between the lower lobe of that curve and the x axis, this very much impacts how I would want to approach the problem.

The scheme I showed you how to implement will be pretty robust to approximate the area contained within the closed curve. That is because it is insensitive to high frequency oscillations in the estimated curve, because those oscillations have high frequency, but small amplitude.

[~,~,~,~,~,~,area] = FSfit(x,y,10,0)
area =
          13979.5620217219
          
>> [~,~,~,~,~,~,area] = FSfit(x,y,20,0)
area =
          13726.1735377685
          
>> [~,~,~,~,~,~,area] = FSfit(x,y,30,0)
area =
          13708.1851809597
          
>> [~,~,~,~,~,~,area] = FSfit(x,y,100,0)
area =
          13707.1667535602

Given the noise, the prediction seems to be fairly stable as long as you use 20 terms in the model, or so.

But to find the area under that curve will be more difficult if we are working with a model of the form I built. So IF my goal was to compute the area blow the curve and above the x axis, I would do it like this:

xmin = min(x);
xmax = max(x);
xadj = [x;xmin;xmax];
yadj = [y;0;0];

E1 = convhull(x,y);
E2 = convhull(xadj,yadj);
polyarea(xadj(E2),yadj(E2)) - polyarea(x(E1),y(E1))
ans =
          4811.42632861809

So we approximate the area under the curve as the difference of two areas, thus the area of the region in purple here:

plot(polyshape(xadj(E2),yadj(E2)))
hold on
plot(polyshape(x(E1),y(E1)))

[Niklas Kurz]

Oh.My.Gosh, you're One Of The Best! I'm like 100000% satisfied now. Excuse me for my randomness, sometimes even I don't know what I'm searching for.

[tarek hussein]

Hello Dear John,

My question is: how i can fit an ellipse from random noisy ring(random data set) as the figure bellow?

i have a polar plot (Nyquist plot) for measuring point on a rotating shaft in X direction( plotting the real and imaginary parts of complex vector X) and i want to find an appropriate curve(circle or ellipse) for that polar curve since i severely need to now the diameter of that ellipse.

I had used the following code to fit:

t=0:0.01:2*pi;

x= a(:,1) % the real part of complex vector which dimension is ( 206310 x1)

y=b(:,1) % the imaginary part of the previous vector

a=pi/4;

z=[cos(a) -sin(a);sin(a) cos(a)];

m=[transpose(x);transpose(y)];

k=z*m;

x=k(1,:);

y=k(2,:);

scatter(x,y);

axis equal;

initialparameter=[12 115 0]; %initialparameter=[x y z]

mx=@(initialparameter)fita(initialparameter,x,y);

outputparameters=fminsearch(mx,initialparameter);

The used function:

function E=fita(initialparameter,xm,ym)

t=linspace(0,2*pi,length(xm));

xao=initialparameter(1)*cos(t);

yao=initialparameter(2)*sin(t);

a=initialparameter(3);

z=[cos(a) -sin(a);sin(a) cos(a)];

m=[xao;yao];

k=z*m;

xao=k(1,:);

yao=k(2,:);

Ex=sum(abs(xao-xm).^2);

Ey=sum(abs(yao-ym).^2);

E=(Ex+Ey)/2;

clf

scatter(xm,ym);

hold on;

plot(xao,yao,'k');

hold off;

axis equal;

drawnow;

end

Note: the code is ok but i have no result appears. the message i have when i run the code is the following:

maximum number of function evaluations has been exceeded - increase MaxfunEvals option. current function value:NaN

Kindly ask you to find a solution of my problem.

Regards

[Image Analyst]

@tarek hussein see the attached paper.

Sorry, I don't have MATLAB implementations for any of the shapes.

[tarek hussein]

thanks for your response