Area of a implicit curve
10 views (last 30 days)
Show older comments
Let's say that I have a multiply connected curve defined implicitly. For example:
f=@(x,y) sin(x).*sin(y)-0.5;
fimplicit(f,[-10,10,-10,10]);
How can we find the area of from the rectangle ?
Here it is explained how can we fill this region: https://www.mathworks.com/matlabcentral/answers/805046-how-can-we-fill-a-region-defined-by-a-implicit-function?s_tid=srchtitle
So, if we can calculate the area of a filled figure than it's all done. There is a function for that kind of evaluations?
P.S. I wrote myself a program which detects where the changes from a connected component to another is made by fimplicit data points and I made a loop that sums all the areas of the polygons created using the [k,v]=boundary(Points,1) command.
I have the areas precision up to 4 decimal places, but I cannot solve the problems that appear on the corners of the rectangle , because boundary don't have them included in the points obtained from fimplicit. To be more precis, if I have the area that my code gives is not the area of a quarter of a circle, but the area of that sector minus the area of the triangle determined by it.
Also I have to mention that the integral method fails (it has too little precision). I mean that we can use the formula
integral2(@(x,y)heaviside(f(x,y)),-10,10,-10,10)
0 Comments
Accepted Answer
Matt J
on 17 Apr 2021
Edited: Matt J
on 18 Apr 2021
Perhaps as follows?
plotRange=[-1,+1, -1,+1]*8;
plotCorners=[-1 1;1 1;1 -1; -1 -1]*8;
f=@(x,y) (sin(x).*sin(y)-0.5);
fp=fimplicit(f,2*plotRange,'MeshDensity',3000);
XY=[fp.XData;fp.YData].'; close
shpRegions=polyshape(XY);
shpBounds=polyshape(plotCorners);
shpRegions=intersect(shpBounds,shpRegions);
Area=area(shpRegions)
plot(shpRegions,'FaceColor','r')
2 Comments
Matt J
on 18 Apr 2021
Here is a modification to detect the complementary area
plotRange=[-1,+1, -1,+1]*8;
plotCorners=[-1 1;1 1;1 -1; -1 -1]*8;
f=@(x,y) (sin(x).*sin(y)-0.5);
fp=fimplicit(f,2*plotRange,'Mesh',4000);
XY=[fp.XData;fp.YData].'; close
shpRegions=polyshape(XY);
shpBounds=polyshape(plotCorners);
shpRegions=intersect(shpBounds,shpRegions);
[x,y]=shpRegions.centroid(1:numboundaries(shpRegions));
s=f(x,y); %test the sign of the centroids
if any( s>0 ) %take complementary region
shpRegions=subtract(shpBounds,shpRegions);
end
Area=area(shpRegions)
plot(shpRegions,'FaceColor','r')
More Answers (1)
Matt J
on 18 Apr 2021
Edited: Matt J
on 18 Apr 2021
Here is a very similar strategy based on alphaShape instead of polyshape. It is slower, but IMO handles the bounding rectangle as well as the distinction between and more gracefully.
plotRange=[-1,+1, -1,+1]*8;
f=@(x,y) (sin(x).*sin(y)-0.5);
fp=fimplicit(f,plotRange,'MeshDensity',3000);
XY=unique([fp.XData;fp.YData].','rows','stable'); close
XY( any( isnan(XY),2) , : )=[];
[Xb,Yb]=ndgrid(linspace(plotRange(1), plotRange(end),500));
supp=f(Xb,Yb)<=0;
shp=alphaShape([XY;Xb(supp), Yb(supp)],0.1);
Area=area(shp)
plot(shp,'FaceColor','r','EdgeColor','none')
2 Comments
See Also
Categories
Find more on Elementary Polygons in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!