http://www.mathworks.com/matlabcentral/newsreader/view_thread/245894
MATLAB Central Newsreader  how to calculate area of interpolated surface
Feed for thread: how to calculate area of interpolated surface
enus
©19942015 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Wed, 04 Mar 2009 14:14:02 +0000
how to calculate area of interpolated surface
http://www.mathworks.com/matlabcentral/newsreader/view_thread/245894#632445
Tomas
I want to compute the area of a interpolated surface. I used the "interp2" method so I got no function relation between the (x,y) and z points. <br>
<br>
Does any of you have a suggestion how to do this?<br>
<br>
Thanks.

Wed, 04 Mar 2009 15:29:01 +0000
Re: how to calculate area of interpolated surface
http://www.mathworks.com/matlabcentral/newsreader/view_thread/245894#632469
John D'Errico
"Tomas " <tomasnic@gmail.com> wrote in message <gom2b9$m0n$1@fred.mathworks.com>...<br>
> I want to compute the area of a interpolated surface. I used the "interp2" method so I got no function relation between the (x,y) and z points. <br>
> <br>
> Does any of you have a suggestion how to do this?<br>
> <br>
> Thanks. <br>
<br>
Since you used interp2, I'll assume that you<br>
have a lattice of points on a regular grid, and<br>
therefore a function evaluated at the lattice of<br>
points f(x,y).<br>
<br>
Now, to compute the area of that surface, there<br>
are several simple solutions. The first is to use<br>
a rectangle rule. So assume we have a function<br>
F as an nxm array, evaluated on the nxm lattice<br>
of points. Assume the spacing in x and y is known<br>
and fixed, at dx and dy respectively.<br>
<br>
Assume the surface is defined as a series of<br>
rectangles, at the height of one of the corners<br>
of each square in the lattice. Then the surface<br>
area is given as<br>
<br>
sarea1 = (n1)*(m1)*dx*dy<br>
<br>
This is the rectangle rule integration. As you can<br>
see, it misses one row and column entirely. It<br>
very much depends on how you will define the<br>
surface though. If instead, you wish to define<br>
the surface as a series of rectangles, each<br>
centered exactly on each lattice point, then the<br>
surface area is simply<br>
<br>
sarea2 = n*m*dx*dy;<br>
<br>
dx = 3;<br>
dy = 2;<br>
n = 10;<br>
m = 20;<br>
F = rand(n,m);<br>
surf(F)<br>
<br>
Most likely, both of these approximations are<br>
not adequate for your purposes. So a better<br>
solution might be to form a set of triangles.<br>
Triangulate the lattice, breaking each square<br>
into a pair of triangles. Find the area of each<br>
triangle, then sum over the areas. It is doable,<br>
and perhaps before I get done writing I'll do<br>
exactly that.<br>
<br>
sarea1 = (n1)*(m1)*dx*dy<br>
sarea1 =<br>
1026<br>
<br>
sarea2 = n*m*dx*dy<br>
sarea2 =<br>
1200<br>
<br>
(Ok, I did it just now, but the method I actually<br>
used was not easy to explain, and would just<br>
have confused things.) Here is the area of the<br>
triangles. As you should expect, it is in the<br>
middle of the others.<br>
<br>
satri<br>
satri =<br>
1170.7<br>
<br>
Here is another estimate, formed by approximating<br>
the surface area by integrating the length of the<br>
piecewise linear line segment along one of the<br>
two dimensions.<br>
<br>
trapz(sum(dx*sqrt(1 + diff(F,[],1).^2),1),2)*dy<br>
ans =<br>
1099<br>
<br>
I'll choose to do the same in the other dimension<br>
too, just to see if it made a difference.<br>
<br>
trapz(sum(dy*sqrt(1 + diff(F,[],2).^2),2),1)*dx<br>
ans =<br>
1101.9<br>
<br>
Happily the difference is slight. The point is,<br>
your solution must depend upon how you define<br>
the surface. There are other, higher order<br>
approximations one might choose for all of this.<br>
<br>
I'll just stop here.<br>
<br>
John

Thu, 05 Mar 2009 16:00:19 +0000
Re: how to calculate area of interpolated surface
http://www.mathworks.com/matlabcentral/newsreader/view_thread/245894#632780
Tomas
Hi John, <br>
<br>
thanks for your answer, but the problem is that I got no function z = F(x,y) at the lattice points. <br>
<br>
Maybe it will help if I describe what I want to do. I want to compute the area of a small region of a topographical map. Maybe I can fit some kind of polynomial function to points on the contour lines using interpolation or a overestimation method. When I got a function I will try the method you described.<br>
<br>
Tomas

Thu, 05 Mar 2009 16:25:04 +0000
Re: how to calculate area of interpolated surface
http://www.mathworks.com/matlabcentral/newsreader/view_thread/245894#632789
John D'Errico
"Tomas " <tomasnic@gmail.com> wrote in message <goosuj$nr0$1@fred.mathworks.com>...<br>
> Hi John, <br>
> <br>
> thanks for your answer, but the problem is that I got no function z = F(x,y) at the lattice points. <br>
> <br>
> Maybe it will help if I describe what I want to do. I want to compute the area of a small region of a topographical map. Maybe I can fit some kind of polynomial function to points on the contour lines using interpolation or a overestimation method. When I got a function I will try the method you described.<br>
> <br>
> Tomas<br>
<br>
Then what do you have? All that I assumed is that<br>
you have a list of values, defined at each point on<br>
a regular lattice. If you used interp2, then you MUST<br>
have that much.<br>
<br>
Do you have points only on the contour lines? If<br>
so, then how do you claim to have used interp2?<br>
This would be impossible.<br>
<br>
So tell me what exactly you do have? Be precise<br>
when you ask a question, else all you do is waste<br>
the time of someone who might try to answer the<br>
wrong question (as apparently did I.)<br>
<br>
John

Thu, 05 Mar 2009 19:02:02 +0000
Re: how to calculate area of interpolated surface
http://www.mathworks.com/matlabcentral/newsreader/view_thread/245894#632821
Roger Stafford
"Tomas " <tomasnic@gmail.com> wrote in message <goosuj$nr0$1@fred.mathworks.com>...<br>
> .......<br>
> I want to compute the area of a small region of a topographical map. <br>
> .......<br>
<br>
Tomas, I am guessing you want to find the approximate surface area of a region on your map, as opposed to the area of that region as projected onto a horizontal plane. That is, you want the "slant" area which takes into account the gradients present in the region. And you wish to use the output of 'interp2' for this purpose.<br>
<br>
It should be pointed out that in the case of mountainous terrains such a surface figure would be highly dependent on the fineness of resolution of your data. Down at a fine resolution the surface area may be vastly greater than at a course resolution because of various surface irregularities. (As a onetime backpacker I can testify to the truth of that statement!)<br>
<br>
In calculus for an infinitesimal rectangular region dx by dy, its threedimensional surface area is given by<br>
<br>
sqrt(1 + (dz/dx)^2+(dz/dy)^2)*dx*dy<br>
<br>
You can therefore use finite differences for approximating this expression for each small rectangular region in the mesh of points XI,YI,ZI to get a first order approximation to the desired surface area. Assuming your data has a reasonably high resolution, that should be adequate for your purposes.<br>
<br>
Roger Stafford

Thu, 05 Mar 2009 19:41:01 +0000
Re: how to calculate area of interpolated surface
http://www.mathworks.com/matlabcentral/newsreader/view_thread/245894#632828
Roger Stafford
"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gom6nt$8m1$1@fred.mathworks.com>...<br>
> .......<br>
> trapz(sum(dy*sqrt(1 + diff(F,[],2).^2),2),1)*dx<br>
> ans = 1101.9<br>
> .......<br>
> John<br>
<br>
Note to John: I see you have also done a surface calculation above but I think that gets an answer that is too small. Shouldn't it be changed to something like this:<br>
<br>
trapz(sum(dy*sqrt(1 + diff(F,[],1).^2 + diff(F,[],2).^2),2),1)*dx<br>
<br>
to get an accurate multiplication by the true secant of the angle between the vertical and the surface normal? Or is my head not working properly this morning?<br>
<br>
Roger Stafford

Thu, 05 Mar 2009 21:36:01 +0000
Re: how to calculate area of interpolated surface
http://www.mathworks.com/matlabcentral/newsreader/view_thread/245894#632858
Tomas
"John D'Errico" <woodchips@rochester.rr.com> wrote in message <br>
> Then what do you have? All that I assumed is that<br>
> you have a list of values, defined at each point on<br>
> a regular lattice. If you used interp2, then you MUST<br>
> have that much.<br>
> <br>
> Do you have points only on the contour lines? If<br>
> so, then how do you claim to have used interp2?<br>
> This would be impossible.<br>
> <br>
> So tell me what exactly you do have? Be precise<br>
> when you ask a question, else all you do is waste<br>
> the time of someone who might try to answer the<br>
> wrong question (as apparently did I.)<br>
> <br>
> John<br>
<br>
I?m very sorry wasting your time, John. I realize now that I can?t use interp2 because i don?t have function values for an equidistant grid. I promise to think more about my questions before posting later.<br>
<br>
Roger, thanks for your answer. I now have an idea what to do.<br>
<br>
Tomas

Thu, 05 Mar 2009 22:23:02 +0000
Re: how to calculate area of interpolated surface
http://www.mathworks.com/matlabcentral/newsreader/view_thread/245894#632871
John D'Errico
"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gop9sd$3ib$1@fred.mathworks.com>...<br>
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <gom6nt$8m1$1@fred.mathworks.com>...<br>
> > .......<br>
> > trapz(sum(dy*sqrt(1 + diff(F,[],2).^2),2),1)*dx<br>
> > ans = 1101.9<br>
> > .......<br>
> > John<br>
> <br>
> Note to John: I see you have also done a surface calculation above but I think that gets an answer that is too small. Shouldn't it be changed to something like this:<br>
> <br>
> trapz(sum(dy*sqrt(1 + diff(F,[],1).^2 + diff(F,[],2).^2),2),1)*dx<br>
<br>
Not quite, though I was sloppy. Best is probably this:<br>
<br>
[X,Y] = meshgrid(0:.1:1,0:.2:3);<br>
dx = 0.1;<br>
dy = 0.2;<br>
[Fx,Fy] = gradient(F,dx,dy);<br>
trapz(trapz(sqrt(1+Fx.^2 + Fy.^2)))*dx*dy<br>
ans =<br>
46.689<br>
<br>
Which correctly integrates the area, using gradient<br>
to compute the partials.<br>
<br>
This generates a surface area that is quite consistent<br>
with an area computed from the sum of the triangulated<br>
area as 46.811.<br>
<br>
I don't know how to solve the problem as posed by the<br>
OP, since we still don't know what he has.<br>
<br>
John

Sat, 07 Mar 2009 08:38:02 +0000
Re: how to calculate area of interpolated surface
http://www.mathworks.com/matlabcentral/newsreader/view_thread/245894#633195
Tomas
"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gopjc6$4i0$1@fred.mathworks.com>...<br>
> I don't know how to solve the problem as posed by the<br>
> OP, since we still don't know what he has.<br>
<br>
I got a topographic map with contour lines. I want to approximate the area of a region of the map based on the height of the contour lines.<br>
<br>
Tomas

Sat, 07 Mar 2009 09:37:01 +0000
Re: how to calculate area of interpolated surface
http://www.mathworks.com/matlabcentral/newsreader/view_thread/245894#633198
Roger Stafford
"Tomas " <tomasnic@gmail.com> wrote in message <gotbpa$prr$1@fred.mathworks.com>...<br>
> I got a topographic map with contour lines. I want to approximate the area of a region of the map based on the height of the contour lines.<br>
> .......<br>
<br>
I'm just speculating here, but suppose you have a way of obtaining the varying "orthogonal" distance between adjacent contour lines. If h is the fixed height difference represented between the contours and d the orthogonal horizontal distance represented between them, then the slant distance across between the contours is just sqrt(d^2+h^2). Integrating this varying slant distance along the space between the contours gives the surface area it represents, as opposed to the horizontally projected area that appears on your map between the contours. Repeating this for all spacings between contours for a region would give you its total surface area.<br>
<br>
Of course this is only approximate, but then being restricted to discrete contour lines is necessarily an inexact procedure. However I am not sure just how you would go about determining this variable quantity d above starting with just a contour map. If the contours are close together they would ordinarily be nearly parallel but widelyspaced contours can sometimes be decidedly nonparallel. Fortunately this latter is just the case where the needed correction over projected area is the smallest.<br>
<br>
Well this is just a wild idea for you to play with.<br>
<br>
Roger Stafford

Sat, 07 Mar 2009 09:48:01 +0000
Re: how to calculate area of interpolated surface
http://www.mathworks.com/matlabcentral/newsreader/view_thread/245894#633199
Bruno Luong
"Tomas " <tomasnic@gmail.com> wrote in message <gotbpa$prr$1@fred.mathworks.com>...<br>
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <gopjc6$4i0$1@fred.mathworks.com>...<br>
> > I don't know how to solve the problem as posed by the<br>
> > OP, since we still don't know what he has.<br>
> <br>
> I got a topographic map with contour lines. I want to approximate the area of a region of the map based on the height of the contour lines.<br>
> <br>
> Tomas<br>
<br>
As Roger pointed out, we need to be able to evalute the gradient of the topographic map at any point in order to compute the area. Therefore I would suggest to approximate your surface map by the spline surface  it is twice derivable  so it is nice to integrate numerically its gradient. That might be violate the property of the map, because land might have vertical fall (infinity gradient), up to you...<br>
<br>
To find the spline surface, you could create a regular grid. Next find find out what to fill on the grid in order to have the contour line data. This is a linear problem : you can solve it by building the matrix (using interp2) and using "\" operator.<br>
<br>
Now there is a little bit of work to extract the gradient at a given point. You can either do your own calculation from the spline data values. You might want to know Matlab use a notaknot condition for interpolation to calculate for points at the outer ring. You could also use 1D spline functions applied to 1D grid crossing the point of interest (that gives you a polynomial coefficient on all braket that you can take the derivative). I believe the spline toolbox have some nice feature to do this kind of calculation.<br>
<br>
The last step: The calculation of the area by using function quad2d. You have to program a function that returns 0 outside the region, and in infinitisimal surface area (see Roger's above formula) inside the region  now you have everything to calculate it. You could also do a brute force sum of infinitisimal surface areas with point that distributed regularly inside the region. You might need a tool to tell when a point is inside the polygonal contour line. I believe there are several such functions on FEX. Be aware about the regions with hole inside it that need some care in your implementation.<br>
<br>
Good luck<br>
<br>
Bruno