Area of a implicit curve
9 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Maxim Bogdan
el 17 de Abr. de 2021
Comentada: Maxim Bogdan
el 18 de Abr. de 2021
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 comentarios
Respuesta aceptada
Matt J
el 17 de Abr. de 2021
Editada: Matt J
el 18 de Abr. de 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 comentarios
Matt J
el 18 de Abr. de 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')
Más respuestas (1)
Matt J
el 18 de Abr. de 2021
Editada: Matt J
el 18 de Abr. de 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 comentarios
Ver también
Categorías
Más información sobre Elementary Polygons en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!