How to apply a fit on this shape?

Question

Niklas Kurz 2020-7-15

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/565277-how-to-apply-a-fit-on-this-shape

移动：John D'Errico 2023-2-20

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)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

John D'Errico 2020-7-15

2
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/565277-how-to-apply-a-fit-on-this-shape#answer_466013

编辑：John D'Errico 2020-7-15

在 MATLAB Online 中打开

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.

9 个评论
显示 7更早的评论隐藏 7更早的评论

John D'Errico 2020-7-15

编辑：John D'Errico 2020-7-19

在 MATLAB Online 中打开

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.

John D'Errico 2020-7-19

在 MATLAB Online 中打开

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

tarek hussein 2023-2-18

移动：John D'Errico 2023-2-20

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 2023-2-19

移动：John D'Errico 2023-2-20

Least Squares orthogonal distances fitting of circle, sphere, ellipse, hyperbola, and parabola.pdf

@tarek hussein see the attached paper.

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

tarek hussein 2023-2-20

移动：John D'Errico 2023-2-20

thanks for your response

请先登录，再进行评论。

How to apply a fit on this shape?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

9 个评论
显示 7更早的评论隐藏 7更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

How to apply a fit on this shape?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

9 个评论 显示 7更早的评论隐藏 7更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

9 个评论
显示 7更早的评论隐藏 7更早的评论