Compute that portion of a surface above a threshold

Hi. I have a variable Z define over a grid of points X,Y. I would like to compute the area of my domain where Z>threshold.

Right now I tried to compute it as

A = trapz(X_vector, trapz(Y_vector, Z> threshold));

This gives acceptable results, however this approach set the border between the "in" and the "out" areas always in the middle of the last point outside and the first point inside.

To better explain the problem, consider the 1-dimensional case where

 x = 1:10;
 y = ones(1,10);
 y(5) = 2;
 threshold = 1.99;
 plot(x,y)
 hold all
 plot(xlim,[1 1]*threshold ,'--k')

Analytically, the portion of the curve above the threshold is equal to 0.2 x-units. However if we plot

 plot(x,y>threshold)

and we compute

 L = trapz(x,y>threshold)
 L =

     1

This is due to the fact that this approach ignores that y(5) is way closer to the threshold than y(4) and y(6).

On the 2D case You can visualize the "correct" solution using

contour(X,Y,Z,[1 1]*threshold)

The problem is that if the contour line is open the computation of the area becomes tricky.

Any suggestions?

Thank you

Answer 1

[Ahmet Cecen]

The problem here is that you appear to be committing a fundamental fallacy in how discrete data work. In your example case, you are plotting this:

and claiming the "real" line segment under the threshold is 0.2 units long. However what you are actually inputting as data looks more like this:

in which case you can clearly see why MATLAB thinks the line segment is 1 units long. You are treating the graphical consequence of MATLAB's plotting by connecting the dots with lines as if it is the actual nature of the data.

I am assuming you are here because finding the analytical definition and integrating for your actual problem is a headache or not possible. I would suggest one of the following:

1) If you have a rough idea about the underlying nature of the surface, you can try to fit an analytical form to it, then find the area above the threshold via integration. For this you should be able to get away with `cftool` function.

2) If your data is just plain ugly, you can use an interpolant to estimate the area. Assuming your data is gridded like your example, use `griddedinterpolant` (`scatteredinterpolant` for ungridded) with an interpolation scheme appropriate for your case. You have treated the above example as if the data points are connected by straight lines, so you can choose `linear` (default) for that. Now you sample your data much denser, for the above example -> `Denser = F(4:0.01:6)`. F is the interpolant, and the step 0.01 you can change based on your accuracy needs and computation power. Now you can threshold `Denser` and use trapz, simple pixel count, a triangulation method or something like marching squares/cubes to get the area.

[Luca Amerio]

I see what you are saying. I was indeed imply to use a linear (or better) interpolation for the data. Isn't this what contour does?

The "use finer grid" idea already came in my mind, but this will increase the computational power of the operation. Since I need to repeat it something between 1'000 and 100'000 times, the computational cost matters in my case.

The fact is that the contour algorithm already computes the correct border for the area I want to compute, but it does not enforce the area to be closed. Therefore methods as polyarea fails when applied to the contour matrix/matrixes.

[Ahmet Cecen]

Hmm... If the nature of your data permits this without ambiguity, you can just pad the data with a layer of 0's so that your Contour lines are forced to be closed.

The polyarea over the C that Contour returns actually works accurately for me in 2016b even when its not closed. Definitely much easier than the interpolant, although I suspect `contour` uses it in the background.

I am pondering what else can be done otherwise on the side.