[R] volume of an irregular grid

[(Ted Harding)]

On 03-Dec-03 Karim Elsawy wrote:
> I have a 3d irregular grid of a surface (closed surface)
> I would like to calculate the volume enclosed inside this surface
> can this be done in R
> any help is very much appreciated
> best regards
> karim

Hi Karim,

You should be able to create your own function for this fairly
straightforwardly, on the following lines.

First you will need to pre-process the data of the grid on the
surface. I assume that the grid is equivalent to approximating
the surface as a set of planar facets.

First, determine for each facet the coordinates (X,Y,Z) of a
point (perhaps, conveniently, its centroid) in the plane of the
facet.

Next, determine the area A of each facet.

Next, determine for each facet the direction of the normal to
the plane of the facet which points _out of_ the surface,
as a unit vector (L,M,N).

You can then set up a vector A and two matrices:
  P with columns (X,Y,Z), and
  D with columns (L,M,N).
each with as many rows as there are facets.

Now take an arbitrary point C; but, for best accuracy, it should
be well within the surface and could be at the centroid of the
facets:

  C<-colMeans(A*X)

Now, for each facet, you want the line joining C to P=(X,Y,Z):

  CP<-(X,Y,X)-C

and also as a unit vector:

  CP1<-CP/sqrt(sum(CP*CP))

Next get the cosine of the angle between CP and the normal
D=(L,M,N):

  gamma<-sum(D*CP1)

Finally, sum

  (1/3)*gamma*A*sqrt(sum(CP*CP))

over all the facets.

This will give you the volume contained within the surface
constituted by the facets.

NB: the above outline calculations refer to each facet separately,
i.e. as if taking the matrices row by row. Vectorising it so as to
do the whole thing in one pass, without loops, should be a nice
exercise in the use of 'sweep'

Hoping this helps,
Ted.


--------------------------------------------------------------------
E-Mail: (Ted Harding) <Ted.Harding at nessie.mcc.ac.uk>
Fax-to-email: +44 (0)870 167 1972
Date: 04-Dec-03                                       Time: 22:09:06
------------------------------ XFMail ------------------------------