I want to fit curve on this polar plot around these points,As shown in the image,Can anyone help me please

Please can anyone tell me how to do it?Any way will be fine

[Walter Roberson]

Is the group at the bottom negative r or large theta?

The ones near the origin, are those the same theta for a few points or is the scale just not large enough to show that they are different theta?

[L K]

The bottom are the negative theta values

and at the origin they are different values, I have attached another image with a larger scale(refer img5.1) and img 5.2 is with larger scale for negative theta

[Walter Roberson]

Can you post the actual r and theta so we can experiment?

[L K]

Yes ,here it is...it is like at that particular radius the theta is plotted

theta=88,89,90,91,92,94,96,94,90,89,-100,-102,-104,-105,-104,-102,-101,-100

radius=5,7,11,17,26,39,46,44,32,3,0,18,34,32,33,29,28,20

i want an envelope kind of thing around those points both up(for positive values) and down(for negative).. (by envelope i only mean a curve looking like an envelope)

I am attaching a figure.. the dots represent some angles at some radius.... so i want that envelope around those angles by using curve fitting..

[Image Analyst]

Do you mean like a pie-shaped sector containing the "bounding box" of all the points? Or do you mean the convex hull?

[L K]

I mean the envelope around those black dots....not the circle at origin

[Walter Roberson]

Is it correct that you want to find an equation in theta that generates approximately radius, together with some bounding values creating a curved box near the projected line, such that the curved box has to contain some specified fraction of the original points? Sort of like if you had computed a line and drawn it on thickly and the thick line had to cover most of the points?

If so, then what portion of the original points need to be covered? 1 standard deviation for example?

[Image Analyst]

What black dots? Do you mean the blue stars? And the envelope would be like if you lasooed the blue stars by hand tracing around them? Or used an active contour or alpha shape?

Answer 1

[Image Analyst]

Try this:

theta=[88,89,90,91,92,94,96,94,90,89,-100,-102,-104,-105,-104,-102,-101,-100];
radius=[5,7,11,17,26,39,46,44,32,3,0,18,34,32,33,29,28,20];
polarplot(theta, radius, 'b*');
hold on;
[x, y] = pol2cart(theta, radius);
k = convhull(x, y);
xch = x(k);
ych = y(k);
[thetaCH, rhoCH] = cart2pol(xch, ych);
polarplot(thetaCH, rhoCH, 'ro-');

[L K]

Thank you very much !!!! Indeed a great work...!!! Helped me really :) i tried to figure out this lot many times, you made my job easier

[L K]

can this be done automatically??? as in if i am given any set of data points...will it work in a similar way?

[Image Analyst]

What's your definition of "automatically"? What about my code is not "automatic". There was no manual interaction needed like asking the user a question or asking him to draw anything. Please give the difference between automatic and manual.

[L K]

I was just trying it on some example data points...

here is the code:

VarName1 is just data points from 1 to 10

x=VarName1;
y = abs(sqrt(x));
[xx,yy] = pol2cart(x,y);
k = convhull(xx, yy);
xch = xx(k);
ych = yy(k);
[xCH, yCH] = cart2pol(xch, ych);
plot(xx(k), yy(k), 'ro-',x,y,'b*');

convex hull is like drawing lines on the edge points going outwards unlike concave hull... from what i understood

the output i got is :

[L K]

OK I GOT MY FAULT, i had converted it from poltocart,it was not needed as they were not in degrees, my bad!!

[L K]

This was convex hull which is shown...can it be done by using concave hull????

[Walter Roberson]

See https://www.mathworks.com/help/matlab/ref/boundary.html

Requires R2014b or later.

For releases before that look in the File Exchange for "alpha shapes"

[L K]

if i use k = boundary(x, y);

and keep the rest of the code same,I am not getting the output, it gives me an error saying

Error using alphaShape The input points must be 2D or 3D coordinates in numpoints-by-ndim format.

I think it is because i am plotting the points on polar plot ,and then drawing the boundary.

So what change in the code is necessary ?

[Image Analyst]

Is your point set different than when I took them to create my answer? What is your current point set now?

[L K]

What you showed me is convex hull...

I was trying to use the same data points,for concave hull,by using k = boundary(x, y)

but i am not getting concave hull...please tell me what changes are needed

[Walter Roberson]

Please post your current code.

At the time of the call, are your x and y scalars or are they vectors or are they arrays?

[L K]

theta=[88,89,90,91,92,94,96,94,90,89,-100,-102,-104,-105,-104,-102,-101,-100];
radius=[5,7,11,17,26,39,46,44,32,3,0,18,34,32,33,29,28,20];
polar(theta, radius, 'b*');
hold on;
[x, y] = pol2cart(theta, radius);
k = boundary(x, y);
xch = x(k);
ych = y(k);
[thetaCH, rhoCH] = cart2pol(xch, ych);
polar(thetaCH, rhoCH, 'ro-');

here is the code

[Walter Roberson]

Remember that your theta is not in radians.

theta=[88,89,90,91,92,94,96,94,90,89,-100,-102,-104,-105,-104,-102,-101,-100];
radius=[5,7,11,17,26,39,46,44,32,3,0,18,34,32,33,29,28,20];
theta_rad = theta * pi/180;
polar(theta_rad, radius, 'b*');
hold on;
[x, y] = pol2cart(theta_rad, radius);
k = boundary(x(:), y(:));     %needs column vectors
xch = x(k);
ych = y(k);
[thetaCH, rhoCH] = cart2pol(xch, ych);
polar(thetaCH, rhoCH, 'ro-');

[L K]

thanx :)

[L K]

Is there any way to give a tolerance say of about 10% to 30% for example around this plot ? so that i can say the points which lie within this tolerance is of so-n-so and which lie outside it is so-n so?

[L K]

Hi Image analyst, Can you please help me in how to give a tolerance of 10 or 20% around this plot?

So that i can tell,values which lie in this tolerance belongs to some specific data and outside it to some other data?

[Image Analyst]

How about you compute the centroid of the region, and then the radius of all the vertex points, and then just send them outwards 10 or 20%? That's what I'd do.

[Walter Roberson]

The poster asked a new Question about this; I posted a solution there. It is a small nuisance because it is polar work.

[L K]

The boundary function which i have used to fit a curve,follows what algorithm? or any formula? i needed it badly in order to understand how the boundary function actually goes about , please if any one of you know,,could you help me out?

[Walter Roberson]

The boundary() function uses Alpha Shapes

However, it is my opinion that the result you are coming up with has no useful meaning.

[L K]

which result you mean ?

[Walter Roberson]

The enclosing boundary that you are coming up with has no useful meaning that I have been able to come up with. I see no reason to use linear boundaries on a curved surface.

[L K]

But is there any way in which i can join the points on the curved surface ?Basically i need to fit a curve around those points... so i thought this was quite closely satisfying the need... but how can i form a curved boundary if not linear boundaries?

[Image Analyst]

Why does it matter? What are you going to do with a boundary, either curved or linear, if you had it? Maybe you think you need it for some reason but actually you don't. If you tell us what you're going to do with the information, maybe we can say whether that is the best approach or not.

[Walter Roberson]

What mathematical properties does the boundary need to have?

[L K]

HI Image analyst, i just have some data points..which i will plot on polar graph and fit boundary around those points...so the next time i run the same thing and if my data points tend to fall outside i will reject it or if within i accept it...this is just overview..... do u have some opinion on this ?Please feel free to share the same ...or if you have some better idea...

[L K]

About mathematical boundaries i am not sure,basically what i know is i need to form an curve like around the points...like i mentioned above ..but after the curve is formed ,,is there a way to know its model ?

Answer 2

[Walter Roberson]

If you have the curve fitting toolbox then use cftool. Import your theta and r and play with some models. Or use pol2cart first and import the resulting x and y and fit those.

[Walter Roberson]

My work so far suggests that a reasonable fit is:

R = -67506.3123723029 - 425.950909132369*Theta + 25140.4082290502*Theta^2 + 86472.2200845025 * sin(Theta + 1.56195303516666) + 6006.21747802037 * sin(Theta^2 + 4.93132363492489)

Here, Theta is theta * Pi / 180 -- that is, theta converted to radians.

The model I fitted against was

a0 + a1 * t + b1 * sin(t + c1) + a2 * t^2 + b2 * sin(t^2 + c2)

Here, without loss of generality, c1 and c2 can be restricted to 0 to 2*Pi .

My fitting was showing two distinct areas of solution, which was really slowing it down. After some thought I realized that because sin(theta+Pi) = -sin(theta), that if b1 and c1 are solutions then so are -b1 and c1+Pi. Therefore without loss of generality you can restrict b1 and b2 to be non-negative as well.

[Walter Roberson]

Marginally better:

-67507.8055608757 - 425.954107825018*Theta + 25140.966254751*Theta^2 + 86474.0893419859*sin(Theta + 1.56195316282131) + 6006.33822370626*sin(Theta^2 + 4.93132643353749)

The range over which the residue values are "quite close" to the above are

-67508.4024493893 -425.955230769825 25140.7891946898 86473.496119004  1.56195312294972 6006.29986670352 4.93132555131309

to

-67507.3317487232 -425.953060425301 25141.1894099114 86474.8375860513 1.56195321779578 6006.38638459066 4.93132751208825

so there is some free play on a0, a2, b1, and a slight bit on b2. The free play is largest on b1 and a0.

The plot looks like this, where blue is the projected values and red is the actual values.

[Walter Roberson]

Visually you do not do all that badly with

- 8.202817146914601 + 31.00545829304457*Theta + 299.819287460388*sin(Theta + 4.628895715827117)

again Theta = theta * pi/180

The other model has a better fit.

This is the model

a0 + a1 * Theta + b1 * sin(Theta + c1)

If you construct

   syms a0 a1 b1 c1 f2(t)
   f2(t) = a0 + a1*t + b1*sin(t + c1)
   obj2 = sum((f2(theta*Pi/180) - radius).^2)
   CCV2 = matlabFunction(obj2, 'Vars', {[a0, a1, b1 c1]})
   fmincon(CCV2, rand(1,4)*2, [], [], [], [], [-inf -inf 0 0], [inf inf inf 2*pi])

then you will probably get the above result. If not, then try the fmincon again; sometimes the rand() puts it into a worse location

[L K]

yes i am getting the same prob, but when i tried it again the error went..but am still not getting the above figure... is there some statement to be put?

[Walter Roberson]

   syms a0 a1 b1 c1 f2(t) Pi
   Pi = sym('pi');
   f2(t) = a0 + a1*t + b1*sin(t + c1)
   obj2 = sum((f2(theta*Pi/180) - radius).^2)
   CCV2 = matlabFunction(obj2, 'Vars', {[a0, a1, b1 c1]})
   options = struct('MaxFunEvals', 3000);
   fmincon(CCV2, rand(1,4)*2, [], [], [], [], [-inf -inf 0 0], [inf inf inf 2*pi], [], options)

You might need to repeat the fmincon a small number of times to get the solution I posted above.

[L K]

I repeated the fmincon,but still the same error...i am not understanding why am i not getting that result..

[Walter Roberson]

theta = [88,89,90,91,92,94,96,94,90,89,-100,-102,-104,-105,-104,-102,-101,-100];
  radius = [5,7,11,17,26,39,46,44,32,3,0,18,34,32,33,29,28,20];
  syms a0 a1 b1 c1 f2(t) Pi
   Pi = sym('pi');
   f2(t) = a0 + a1*t + b1*sin(t + c1)
   obj2 = sum((f2(theta*Pi/180) - radius).^2)
   CCV2 = matlabFunction(obj2, 'Vars', {[a0, a1, b1 c1]})
   options = struct('MaxFunEvals', 3000);
   fmincon(CCV2, rand(1,4)*2, [], [], [], [], [-inf -inf 0 0], [inf inf inf 2*pi], [], options)
f2(t) =
a0 + a1*t + b1*sin(c1 + t)
obj2 =
(a0 + (pi*a1)/2 + b1*sin(c1 + pi/2) - 11)^2 + (a0 - (5*pi*a1)/9 + b1*sin(c1 - (5*pi)/9))^2 + (a0 + (pi*a1)/2 + b1*sin(c1 + pi/2) - 32)^2 + (a0 - (5*pi*a1)/9 + b1*sin(c1 - (5*pi)/9) - 20)^2 + (a0 - (7*pi*a1)/12 + b1*sin(c1 - (7*pi)/12) - 32)^2 + (a0 + (8*pi*a1)/15 + b1*sin(c1 + (8*pi)/15) - 46)^2 + (a0 - (17*pi*a1)/30 + b1*sin(c1 - (17*pi)/30) - 18)^2 + (a0 - (17*pi*a1)/30 + b1*sin(c1 - (17*pi)/30) - 29)^2 + (a0 + (22*pi*a1)/45 + b1*sin(c1 + (22*pi)/45) - 5)^2 + (a0 + (23*pi*a1)/45 + b1*sin(c1 + (23*pi)/45) - 26)^2 + (a0 - (26*pi*a1)/45 + b1*sin(c1 - (26*pi)/45) - 33)^2 + (a0 - (26*pi*a1)/45 + b1*sin(c1 - (26*pi)/45) - 34)^2 + (a0 + (47*pi*a1)/90 + b1*sin(c1 + (47*pi)/90) - 39)^2 + (a0 + (47*pi*a1)/90 + b1*sin(c1 + (47*pi)/90) - 44)^2 + (a0 + (89*pi*a1)/180 + b1*sin(c1 + (89*pi)/180) - 3)^2 + (a0 + (89*pi*a1)/180 + b1*sin(c1 + (89*pi)/180) - 7)^2 + (a0 + (91*pi*a1)/180 + b1*sin(c1 + (91*pi)/180) - 17)^2 + (a0 - (101*pi*a1)/180 + b1*sin(c1 - (101*pi)/180) - 28)^2
CCV2 =
  function_handle with value:
    @(in1)(in1(:,1)+in1(:,2).*pi.*(1.0./2.0)+in1(:,3).*sin(in1(:,4)+pi.*(1.0./2.0))-1.1e1).^2+(in1(:,1)-in1(:,2).*pi.*(5.0./9.0)+in1(:,3).*sin(in1(:,4)-pi.*(5.0./9.0))).^2+(in1(:,1)+in1(:,2).*pi.*(1.0./2.0)+in1(:,3).*sin(in1(:,4)+pi.*(1.0./2.0))-3.2e1).^2+(in1(:,1)-in1(:,2).*pi.*(5.0./9.0)+in1(:,3).*sin(in1(:,4)-pi.*(5.0./9.0))-2.0e1).^2+(in1(:,1)-in1(:,2).*pi.*(7.0./1.2e1)+in1(:,3).*sin(in1(:,4)-pi.*(7.0./1.2e1))-3.2e1).^2+(in1(:,1)+in1(:,2).*pi.*(8.0./1.5e1)+in1(:,3).*sin(in1(:,4)+pi.*(8.0./1.5e1))-4.6e1).^2+(in1(:,1)-in1(:,2).*pi.*(1.7e1./3.0e1)+in1(:,3).*sin(in1(:,4)-pi.*(1.7e1./3.0e1))-1.8e1).^2+(in1(:,1)-in1(:,2).*pi.*(1.7e1./3.0e1)+in1(:,3).*sin(in1(:,4)-pi.*(1.7e1./3.0e1))-2.9e1).^2+(in1(:,1)+in1(:,2).*pi.*(2.2e1./4.5e1)+in1(:,3).*sin(in1(:,4)+pi.*(2.2e1./4.5e1))-5.0).^2+(in1(:,1)+in1(:,2).*pi.*(2.3e1./4.5e1)+in1(:,3).*sin(in1(:,4)+pi.*(2.3e1./4.5e1))-2.6e1).^2+(in1(:,1)-in1(:,2).*pi.*(2.6e1./4.5e1)+in1(:,3).*sin(in1(:,4)-pi.*(2.6e1./4.5e1))-3.3e1).^2+(in1(:,1)-in1(:,2).*pi.*(2.6e1./4.5e1)+in1(:,3).*sin(in1(:,4)-pi.*(2.6e1./4.5e1))-3.4e1).^2+(in1(:,1)+in1(:,2).*pi.*(4.7e1./9.0e1)+in1(:,3).*sin(in1(:,4)+pi.*(4.7e1./9.0e1))-3.9e1).^2+(in1(:,1)+in1(:,2).*pi.*(4.7e1./9.0e1)+in1(:,3).*sin(in1(:,4)+pi.*(4.7e1./9.0e1))-4.4e1).^2+(in1(:,1)+in1(:,2).*pi.*(8.9e1./1.8e2)+in1(:,3).*sin(in1(:,4)+pi.*(8.9e1./1.8e2))-3.0).^2+(in1(:,1)+in1(:,2).*pi.*(8.9e1./1.8e2)+in1(:,3).*sin(in1(:,4)+pi.*(8.9e1./1.8e2))-7.0).^2+(in1(:,1)+in1(:,2).*pi.*(9.1e1./1.8e2)+in1(:,3).*sin(in1(:,4)+pi.*(9.1e1./1.8e2))-1.7e1).^2+(in1(:,1)-in1(:,2).*pi.*(1.01e2./1.8e2)+in1(:,3).*sin(in1(:,4)-pi.*(1.01e2./1.8e2))-2.8e1).^2

Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3000 (the selected value).

ans =
          22.5101277310399           1.5446077646475          11.4355301636198          4.51809594253506
>> fmincon(CCV2, rand(1,4)*2, [], [], [], [], [-inf -inf 0 0], [inf inf inf 2*pi], [], options)

Local minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the default value of the optimality tolerance,
and constraints are satisfied to within the default value of the constraint tolerance.

<stopping criteria details>

ans =
          23.5818762479881        -0.285742815508204      1.80149433600663e-07         0.424643259717351
>> fmincon(CCV2, rand(1,4)*2, [], [], [], [], [-inf -inf 0 0], [inf inf inf 2*pi], [], options)

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the default value of the step size tolerance and constraints are 
satisfied to within the default value of the constraint tolerance.

<stopping criteria details>

ans =
          23.5818765447779        -0.285742717659878      3.17390351376664e-08          2.06965743037983
>> fmincon(CCV2, rand(1,4)*2, [], [], [], [], [-inf -inf 0 0], [inf inf inf 2*pi], [], options)

Solver stopped prematurely.

fmincon stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 3000 (the selected value).

ans =
          22.4306120895266           1.6798711087781          12.2504648986848           4.5366287077669
>> fmincon(CCV2, rand(1,4)*2, [], [], [], [], [-inf -inf 0 0], [inf inf inf 2*pi], [], options)

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the default value of the step size tolerance and constraints are 
satisfied to within the default value of the constraint tolerance.

<stopping criteria details>

ans =
         -8.20120347875585          31.0148226862162           299.81610540455          4.62883973782363

If you are seeing an error message, please post it.

[L K]

How do i get that plot of the projected value and the actual value...

[Walter Roberson]

coeffs = fmincon(CCV2, rand(1,4)*2, [], [], [], [], [-inf -inf 0 0], [inf inf inf 2*pi], [], options);
F2 = subs(f2, {a0, a1, b0, c1}, {coefs(1), coeffs(2), coeffs(3), coeffs(4)};
sort_theta = sort(theta);
projected_r = double(F2(sort_theta * Pi/180));
polar(sort_theta, projected_r);