You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
I want to fit curve on this polar plot around these points,As shown in the image,Can anyone help me please
10 views (last 30 days)
Show older comments
Please can anyone tell me how to do it?Any way will be fine
8 Comments
Walter Roberson
on 23 Jan 2017
Edited: Walter Roberson
on 23 Jan 2017
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
on 23 Jan 2017
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
L K
on 28 Jan 2017
Edited: L K
on 28 Jan 2017
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
on 28 Jan 2017
Edited: Image Analyst
on 28 Jan 2017
Do you mean like a pie-shaped sector containing the "bounding box" of all the points? Or do you mean the convex hull?
Walter Roberson
on 28 Jan 2017
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
on 29 Jan 2017
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?
Accepted Answer
Image Analyst
on 28 Jan 2017
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-');
27 Comments
L K
on 29 Jan 2017
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
on 29 Jan 2017
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
on 29 Jan 2017
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
on 30 Jan 2017
Edited: Walter Roberson
on 30 Jan 2017
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
on 30 Jan 2017
OK I GOT MY FAULT, i had converted it from poltocart,it was not needed as they were not in degrees, my bad!!
Walter Roberson
on 31 Jan 2017
Requires R2014b or later.
For releases before that look in the File Exchange for "alpha shapes"
L K
on 31 Jan 2017
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
on 31 Jan 2017
Is your point set different than when I took them to create my answer? What is your current point set now?
L K
on 31 Jan 2017
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
on 31 Jan 2017
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
on 1 Feb 2017
Edited: Walter Roberson
on 1 Feb 2017
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
on 1 Feb 2017
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
on 2 Feb 2017
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
on 10 Feb 2017
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
on 10 Feb 2017
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
on 10 Feb 2017
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
on 22 Mar 2017
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
on 22 Mar 2017
However, it is my opinion that the result you are coming up with has no useful meaning.
Walter Roberson
on 22 Mar 2017
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
on 22 Mar 2017
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
on 22 Mar 2017
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.
L K
on 23 Mar 2017
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
on 23 Mar 2017
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 ?
More Answers (1)
Walter Roberson
on 23 Jan 2017
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.
9 Comments
Walter Roberson
on 28 Jan 2017
Edited: Walter Roberson
on 29 Jan 2017
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
on 29 Jan 2017
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
on 29 Jan 2017
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
on 29 Jan 2017
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
on 30 Jan 2017
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
on 31 Jan 2017
I repeated the fmincon,but still the same error...i am not understanding why am i not getting that result..
Walter Roberson
on 31 Jan 2017
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.
Walter Roberson
on 31 Jan 2017
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);
See Also
Categories
Find more on Computational Geometry in Help Center and File Exchange
Tags
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)