Estimate area inside point cloud

Hi, I'm trying to determine the energy consumed during a hysteresis cycle. That should be the area locked inside a point cloud such as the one in the figure, where two 600 element vectors act as border.

I have tried trapz(H,B), what results in a negative solution, so I was wondering if there was a way of maybe defining the point cloud as a curve and integrating it.

All data comes from experimental results, so I don't have a function to use in integral2.

Answer 1

[Star Strider]

‘I have tried trapz(H,B), what results in a negative solution’

The data are still reliable. The negative solution is the result of the direction that trapz integrates the data.

Example —

t = linspace(0, 2*pi, 1000);
x = 2*cos(t);
y = 2*sin(t);
A = trapz(x, y)

A = -12.5663

The area is still , so:

r = 2;
A = pi*r^2

A = 12.5664

And the sign here is irrelevant.

EDIT — (1 May 2021 at 18:00)

Another illustration —

t = linspace(2*pi, 0, 1000);                % Reversed Direction Produces Positive Result
x = 2*cos(t);
y = 2*sin(t);
A = trapz(x, y)

A = 12.5663

.

Answer 2

[John D'Errico]

As long as the curve is a sequence of points defining the perimeter and they are distinct, you can just compute the area using a tool like polyarea. Or you can create a polyshape, and then compute the area of that curve. Be careful though. You cannot apply any such tool (including trapz) to a randomly sequenced point cloud. They MUST be in sequence around the perimeter.

For example:

n = 100;
theta = linspace(0,2*pi,n);
x = cos(theta);
y = sin(theta);
plot(x,y,'o')
axis equal

Now compute the area.

polyarea(x,y)

ans = 3.1395

The area of a circle of unit radius would be pi, so roughly 3.14. This will be a slight underestimate, becuase the polygonal region lies entirely inside the circle. So 3.1395 is decent.

The nice thing is polyarea does not care if the sequnce traverses the circle in a counter-clockwise or clockwise manner. That is, the above sequence went clockwise around the circle. So this next area, which is traversed in a clockwise manner, is also positive.

polyarea(cos(flip(theta)),sin(flip(theta)))

ans = 3.1395

And we could have used a polyshape too.

ps = polyshape(x(1:end-1),y(1:end-1));

I dropped the last point from each of x and y there so polyshape would not get upset, since that point is essentially also the first point.

plot(ps)
axis equal

area(ps)

ans = 3.1395

[Jose Ignacio Menéndez Guillén]

Thank you!