Area of a implicit curve

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)

Answer 1

[Matt J]

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

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

shpBounds=polyshape(plotCorners);

shpRegions=intersect(shpBounds,shpRegions);

Area=area(shpRegions)

Area = 48.6866

plot(shpRegions,'FaceColor','r')

[Maxim Bogdan]

It's a very nice solution! However the program does not make the difference between f or . It returns the same result and the same region. I wonder if it is posible to introduce somehow the condition ... If it is veryfied by a coloured point the the region is that. Else the region is the complementary one, and that should be filled.

[Matt J]

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

Warning: Polyshape has duplicate vertices, intersections, or other inconsistencies that may produce inaccurate or unexpected results. Input data has been modified to create a well-defined polyshape.

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)

Area = 207.3129

plot(shpRegions,'FaceColor','r')

Answer 2

[Matt J]

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)

Area = 207.3125

plot(shp,'FaceColor','r','EdgeColor','none')

[John D'Errico]

Both good answers.

[Maxim Bogdan]

Thanks a lot Matt J! Now the problem of the area and perimeter of a implicitly defined shape is very clear.