Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Integrate area/volume using triscatteredinterp?

Subject: Integrate area/volume using triscatteredinterp?

From: AJP

Date: 1 Mar, 2011 12:38:05

Message: 1 of 11

Say I have some smooth, single valued 2D function that goes like y=f(x) or some equivalent 3D function z=f(x,y), but I only know it at discrete coordinates.

Triscatteredinterp outputs an interpolant F that will provide the value of the function at some location (x) in 2D or (x,y) in 3D. Therefore typing F(x,y) will return an interpolated z value for example.

In theory couldn't this be used in conjunction with an adaptive quadrature algorithm like quad/dblquad to find the area under the curve/surface?

The only problem is quad() / dblquad() need to be passed a function, whereas the iterpolant F, while able to behave like a function, is a TriScatteredInterp class.

A function handle like handle=@F(x,y) cannot be created for this class.

Can anyone think of a trick so that the interpolant F can be passed to quad like a normal function?

Subject: Integrate area/volume using triscatteredinterp?

From: AJP

Date: 1 Mar, 2011 12:55:19

Message: 2 of 11

"AJP" wrote in message <ikipbd$llu$1@fred.mathworks.com>...
> Say I have some smooth, single valued 2D function that goes like y=f(x) or some equivalent 3D function z=f(x,y), but I only know it at discrete coordinates.
===============================================

As an example, say I have:

z=cos(x)+cos(y)

and I want to find the volume under the surface z for -1<x<1 , -1<y<1.

With dblquad I can simply define a function handle and specify the limits to obtain the answer:

I=dblquad(@(x,y)cos(x)+cos(y),-1,1,-1,1)


Now suppose I know z only at discrete (x,y), but I do know that it's smooth and continuous, so things like noise and discontinuities aren't an additional problem.

I can calculate F=TriScatteredInterp(x,y,z), such that now if I put F(x1,y1) I get an interpolated z value.

This means it will work in a similar way to fhandle=(@(x,y)cos(x)+cos(y).


My question is then, is there a way to evaluate something like:

I=dblquad(@F(x,y),-1,1,-1,1)

and obtain a similar result as before with a normal function handle?

Subject: Integrate area/volume using triscatteredinterp?

From: AJP

Date: 1 Mar, 2011 13:19:23

Message: 3 of 11

"AJP" wrote in message <ikiqbn$qqu$1@fred.mathworks.com>...
> "AJP" wrote in message <ikipbd$llu$1@fred.mathworks.com>...
> > Say I have some smooth, single valued 2D function that goes like y=f(x) or some equivalent 3D function z=f(x,y), but I only know it at discrete coordinates.
> ===============================================

Hmm, it seems my problem was because dblquad may at times evaluate the function passed to it with vector x and scalar y.

This is a requirement that the TriScatteredInterp interpolant F cannot satisfy.

e.g.:

F([0,1],1)
??? Error using ==> subsref
The input data has inconsistent size.

However, I have managed to make it work using quad2d.

This is probably preferable anyway since, strictly speaking, one should not use higher-order or adaptive integration algorithms on unknown functions / discrete data.

The solution is therefore as follows:

x = linspace(-1,1);
y = linspace(-1,1);
[X Y]=meshgrid(x,y);
Z = cos(X)+cos(Y);
F = TriScatteredInterp(X(:),Y(:),Z(:));
I1=dblquad(@(x,y)cos(x)+cos(y),-1,1,-1,1)
I2=quad2d(@(a,b)intfunc(F,a,b),-1,1,-1,1)

I1 =
   6.7318

I2 =
   6.7315

Note I have the following in a separate function m file:

function a=intfunc(F,x,y)
a=F(x,y);
end

Therefore an originally scattered dataset can be integrated over a rectangular area.

Now I have to do find a way to handle nonrectangular cases.....!

Subject: Integrate area/volume using triscatteredinterp?

From: Bruno Luong

Date: 1 Mar, 2011 13:33:23

Message: 4 of 11

"AJP" wrote in message <ikipbd$llu$1@fred.mathworks.com>...
> Say I have some smooth, single valued 2D function that goes like y=f(x) or some equivalent 3D function z=f(x,y), but I only know it at discrete coordinates.
>
> Triscatteredinterp outputs an interpolant F that will provide the value of the function at some location (x) in 2D or (x,y) in 3D. Therefore typing F(x,y) will return an interpolated z value for example.
>
> In theory couldn't this be used in conjunction with an adaptive quadrature algorithm like quad/dblquad to find the area under the curve/surface?
>
> The only problem is quad() / dblquad() need to be passed a function, whereas the iterpolant F, while able to behave like a function, is a TriScatteredInterp class.
>
> A function handle like handle=@F(x,y) cannot be created for this class.

Why not? There need no trick here, just follow the doc of dblquad

x = rand(100,1)*4-2;
y = rand(100,1)*4-2;
z = x.*exp(-x.^2-y.^2);

F = TriScatteredInterp(x,y,z);

f = @(x,y) F(x,y+zeros(size(x)))
I = dblquad(f,-1,1,-1,1)

However I must say this is a poor way of integration. The underlying function returned by TriScatteredInterp is known and the integral can be computed much more efficiently than using dblquad.

Bruno

Subject: Integrate area/volume using triscatteredinterp?

From: John D'Errico

Date: 1 Mar, 2011 13:35:19

Message: 5 of 11

"AJP" wrote in message <ikipbd$llu$1@fred.mathworks.com>...
> Say I have some smooth, single valued 2D function that goes like y=f(x) or some equivalent 3D function z=f(x,y), but I only know it at discrete coordinates.
>
> Triscatteredinterp outputs an interpolant F that will provide the value of the function at some location (x) in 2D or (x,y) in 3D. Therefore typing F(x,y) will return an interpolated z value for example.
>
> In theory couldn't this be used in conjunction with an adaptive quadrature algorithm like quad/dblquad to find the area under the curve/surface?
>
> The only problem is quad() / dblquad() need to be passed a function, whereas the iterpolant F, while able to behave like a function, is a TriScatteredInterp class.
>
> A function handle like handle=@F(x,y) cannot be created for this class.
>
> Can anyone think of a trick so that the interpolant F can be passed to quad like a normal function?

There are several issues here.

First of all, triscatteredinterp offers three distinct methods
for interpolation, ALL of which are non-smooth. Were you
to use those methods in an adaptive scheme like quad, it
would have problems, as those schemes are designed for
smooth functions. When they see non-smooth functions,
they tend to throw lots of points down around to resolve
the derivative singularities. This will happen even for the
smoothest of those options, the linear case.

Assuming that you want to integrate the piecewise linear
interpolant under a triscatteredinterp surface, then in the
simple case for 1-d, just use trapz!!!!!! Since trapz gives
the exact solution for a piecewise linear function, you are
done.

In the case of a function of two variables, where the
function is KNOWN only at the vertices of the triangulation,
with an implied linear interpolation across that triangulation,
the solution is still simple enough.

Compute the mean of the values of your function over the
three vertices of each triangle in the triangulation, then
multiply by the area of the corresponding triangle in the
tessellation (in the (x,y) plane). Sum up those results. I
can give you this if you can't figure it out yourself and it
is what you need. In fact, this procedure is generally valid
for the piecewise linear interpolant over a tessellated domain
in any number of dimensions.

If you wanted to compute the volume underneath a general
function over the triangulated domain, this gets slightly
more difficult, but not tremendously so.

For a Gaussian quadrature over that domain, you can use
a function I wrote some time ago to solve the problem in
n-dimensions. tessquad is on the FEX, so find it here:

http://www.mathworks.com/matlabcentral/fileexchange/9991-tessquad

This tool integrates a GENERAL function defined on a
triangulated domain, using a fixed Gaussian rule over each
simplex in that tessellation.

For example, integrate the function z = exp(x+y) over the
unit square. I'll just build a triangulation by hand.

>> format long g
>> fun = @(XY) exp(sum(XY,2));
>> tessquad(fun,2,[0 0;0 1;1 0;1 1],[1 2 3;2 3 4])
ans =
          2.95211130668392

>> tessquad(fun,4,[0 0;0 1;1 0;1 1],[1 2 3;2 3 4])
ans =
          2.95249244121631

>> tessquad(fun,6,[0 0;0 1;1 0;1 1],[1 2 3;2 3 4])
ans =
          2.95249244201256

As you can see, the result has stabilized by the time
I use a rule with reasonable order. Of course, the exact
answer is trivial to write out,

>> (exp(1) - 1)^2
ans =
          2.95249244201256

For the even more general case, where a truly adaptive
rule is desired over a triangulated domain in several
dimensions, I've written code for that too, but I've not
yet gotten my stopping criteria to the point where I am
happy with it.

John

Subject: Integrate area/volume using triscatteredinterp?

From: AJP

Date: 1 Mar, 2011 13:48:04

Message: 6 of 11

Thank you very much for you insightful comments.

I had actually wondered about using a "(mean height) * (area)" type method as John mentions, so I'll give that a go now.

Again thank you to all for your comments.

Subject: Integrate area/volume using triscatteredinterp?

From: AJP

Date: 3 Mar, 2011 15:12:05

Message: 7 of 11

"John D'Errico" <woodchips@rochester.rr.com> wrote in message

> Compute the mean of the values of your function over the
> three vertices of each triangle in the triangulation, then
> multiply by the area of the corresponding triangle in the
> tessellation (in the (x,y) plane). Sum up those results. I
> can give you this if you can't figure it out yourself and it
> is what you need. In fact, this procedure is generally valid
> for the piecewise linear interpolant over a tessellated domain
> in any number of dimensions.

> John

---------------------------------------------------------------------------------------

OK, to recap:

I have a matrix of points, P=[x,y,z] where vectors x and y store my (x,y) points at which I know z (some f(x,y), known only at the given (x,y) positions).

I obtain the Delaunay triangulation using:

tri=delaunay(x,y);

I then compute a vector of triangles areas using the following function with x=P(:,1) and y=P(:,2):

function areas=triareas(x,y,tri)
ab=[x(tri(:,1)),y(tri(:,1))]-[x(tri(:,2)),y(tri(:,2))];
ac=[x(tri(:,1)),y(tri(:,1))]-[x(tri(:,3)),y(tri(:,3))];
areas=abs(ab(:,1).*ac(:,2) - ab(:,2).*ac(:,1))*0.5;

Now I obtain the mean of the three z values at the vertices of each triangle using the following:

function heights=triheights(z,tri)
heights=zeros(length(tri),1);
for n=1:length(tri)
heights(n)=mean(z(tri(n,:)));
end

Then finally the integral is:

int=sum(areas.*heights);

This seems to work fine and answers the question I initially posed. I have two questions:

(1) Is there a way to find the heights (i.e. mean of z over the 3 vertices of each triangle) without using a loop?

(2) Is it my imagination or does the delaunay() function work slower when triangulating N points on a regular grid than when triangulating N scattered points?

Subject: Integrate area/volume using triscatteredinterp?

From: John D'Errico

Date: 3 Mar, 2011 15:42:08

Message: 8 of 11

"AJP" wrote in message <ikob45$cve$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message
>
> > Compute the mean of the values of your function over the
> > three vertices of each triangle in the triangulation, then
> > multiply by the area of the corresponding triangle in the
> > tessellation (in the (x,y) plane). Sum up those results. I
> > can give you this if you can't figure it out yourself and it
> > is what you need. In fact, this procedure is generally valid
> > for the piecewise linear interpolant over a tessellated domain
> > in any number of dimensions.
>
> > John
>
> ---------------------------------------------------------------------------------------
>
> OK, to recap:
>
> I have a matrix of points, P=[x,y,z] where vectors x and y store my (x,y) points at which I know z (some f(x,y), known only at the given (x,y) positions).
>
> I obtain the Delaunay triangulation using:
>
> tri=delaunay(x,y);
>
> I then compute a vector of triangles areas using the following function with x=P(:,1) and y=P(:,2):
>
> function areas=triareas(x,y,tri)
> ab=[x(tri(:,1)),y(tri(:,1))]-[x(tri(:,2)),y(tri(:,2))];
> ac=[x(tri(:,1)),y(tri(:,1))]-[x(tri(:,3)),y(tri(:,3))];
> areas=abs(ab(:,1).*ac(:,2) - ab(:,2).*ac(:,1))*0.5;

This seems right at a glance.


>
> Now I obtain the mean of the three z values at the vertices of each triangle using the following:
>
> function heights=triheights(z,tri)
> heights=zeros(length(tri),1);
> for n=1:length(tri)
> heights(n)=mean(z(tri(n,:)));
> end
>
> Then finally the integral is:
>
> int=sum(areas.*heights);
>

This too.


> This seems to work fine and answers the question I initially posed. I have two questions:
>
> (1) Is there a way to find the heights (i.e. mean of z over the 3 vertices of each triangle) without using a loop?
>

You can do it in a vectorized form if you want. This
should work.

  meanheights = mean(z(tri),2);

> (2) Is it my imagination or does the delaunay() function work slower when triangulating N points on a regular grid than when triangulating N scattered points?

I'd not be surprised. At the same time, triangulating a
regular quadrilateral mesh is trivial to do, and far faster
than delaunay would be anyway. I've got tools for this
on a general lattice in n-dimensions in my simplicial
complex toolbox, but here is a simple 2-d version. For
an n by m array...

n = 2;
m = 3;

[i1,i2] = ndgrid(1:(n-1),1:(m-1));
ind = i1(:) + (i2(:)-1)*n;
tri = [[ind,ind + 1,ind+n+1];[ind,ind+n,ind+n+1]];

>> tri
tri =
     1 2 4
     3 4 6
     1 3 4
     3 5 6

John

Subject: Integrate area/volume using triscatteredinterp?

From: AJP

Date: 3 Mar, 2011 16:17:05

Message: 9 of 11

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <ikocsg$s8l$1@fred.mathworks.com>...

> > (2) Is it my imagination or does the delaunay() function work slower when triangulating N points on a regular grid than when triangulating N scattered points?
>
> I'd not be surprised. At the same time, triangulating a
> regular quadrilateral mesh is trivial to do, and far faster
> than delaunay would be anyway. I've got tools for this
> on a general lattice in n-dimensions in my simplicial
> complex toolbox, but here is a simple 2-d version. For
> an n by m array...
>
> n = 2;
> m = 3;
>
> [i1,i2] = ndgrid(1:(n-1),1:(m-1));
> ind = i1(:) + (i2(:)-1)*n;
> tri = [[ind,ind + 1,ind+n+1];[ind,ind+n,ind+n+1]];
>
> >> tri
> tri =
> 1 2 4
> 3 4 6
> 1 3 4
> 3 5 6
>
> John



Thank you John for your wonderful comments and entertaining my ignorance for so long already.

I shall definitely investigate ways to more quickly mesh my grid. This is currently a problem with delaunay() as I want to work on a fine grid.

Your example is for grid points on a regular, rectangular x,y grid. The problem I'm working on will use such a grid initially, and being able to mesh this quickly is very helpful. Thank you for your help here.

Eventually I would also like to handle a rotated rectangular grid (i.e. an x-y grid rotated about the z axis). The grid remains regular but is now not alligned to the x-y cartesian axes (as in the result rotate(h,[0 0 1],alpha) would have on h=plot(x,y) ).

Have you any tips on going about a fast way of triangulating this?

Perhaps I can use your fast meshing example on the unrotated recangular grid, then rotate the triangulation? My guess would be that this would be easier than trying to mesh an arbitrarily rotated grid.

Anyway, I shall give it a try and see what I can come up with. Thank you for any comments.

Subject: Integrate area/volume using triscatteredinterp?

From: AJP

Date: 3 Mar, 2011 18:18:05

Message: 10 of 11

"AJP" wrote in message <ikoeu1$dmj$1@fred.mathworks.com>...

> Eventually I would also like to handle a rotated rectangular grid (i.e. an x-y grid rotated about the z axis). The grid remains regular but is now not alligned to the x-y cartesian axes (as in the result rotate(h,[0 0 1],alpha) would have on h=plot(x,y) ).
>
> Have you any tips on going about a fast way of triangulating this?
>
> Perhaps I can use your fast meshing example on the unrotated recangular grid, then rotate the triangulation? My guess would be that this would be easier than trying to mesh an arbitrarily rotated grid.
=================================

Having thought about this a little, I think my previous post asked somewhat of a silly question.

Since the triangulation is just a set of indices that reference the points used to create it, then surely I can do any amount of rotation or translation that I like on these - the points' values may change, but they will still be in the same order, hence the triangulation is still valid.

If this is correct then I can indeed perform the triangulation first, then rotate/translate the points and then end up with a triangulated, rotated rectangular grid.

Thank you for the inspiriation!

Subject: Integrate area/volume using triscatteredinterp?

From: John D'Errico

Date: 3 Mar, 2011 19:07:04

Message: 11 of 11

"AJP" wrote in message <ikom0t$m7f$1@fred.mathworks.com>...
> "AJP" wrote in message <ikoeu1$dmj$1@fred.mathworks.com>...
>
> > Eventually I would also like to handle a rotated rectangular grid (i.e. an x-y grid rotated about the z axis). The grid remains regular but is now not alligned to the x-y cartesian axes (as in the result rotate(h,[0 0 1],alpha) would have on h=plot(x,y) ).
> >
> > Have you any tips on going about a fast way of triangulating this?
> >
> > Perhaps I can use your fast meshing example on the unrotated recangular grid, then rotate the triangulation? My guess would be that this would be easier than trying to mesh an arbitrarily rotated grid.
> =================================
>
> Having thought about this a little, I think my previous post asked somewhat of a silly question.
>

Well, perhaps it is best you said that. :-)


> Since the triangulation is just a set of indices that reference the points used to create it, then surely I can do any amount of rotation or translation that I like on these - the points' values may change, but they will still be in the same order, hence the triangulation is still valid.
>
> If this is correct then I can indeed perform the triangulation first, then rotate/translate the points and then end up with a triangulated, rotated rectangular grid.
>

Yep. The triangulation generated is just a set of references.
If you then do any arbitrary rotation or translation on the
point themselves, nothing stops the triangulation from being
valid. Mirror images, any simple transformation are all ok.

In fact, you can even do a bit of nonlinear stuff to the points
while still leaving the transformation a valid one, although at
some point things will fail. (You can identify when that happens.)

John

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us