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?

"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).

"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.

"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.

"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:

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.

> 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.

"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...

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.

> 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.

"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.)

What is a watch list?

You can think of your watch list as threads that you have bookmarked.

You can add tags, authors, threads, and even search results to your watch list. This way you can easily keep track of topics that you're interested in. To view your watch list, click on the "My Newsreader" link.

To add items to your watch list, click the "add to watch list" link at the bottom of any page.

How do I add an item to my watch list?

Search

To add search criteria to your watch list, search for the desired term in the search box. Click on the "Add this search to my watch list" link on the search results page.

You can also add a tag to your watch list by searching for the tag with the directive "tag:tag_name" where tag_name is the name of the tag you would like to watch.

Author

To add an author to your watch list, go to the author's profile page and click on the "Add this author to my watch list" link at the top of the page. You can also add an author to your watch list by going to a thread that the author has posted to and clicking on the "Add this author to my watch list" link. You will be notified whenever the author makes a post.

Thread

To add a thread to your watch list, go to the thread page and click the "Add this thread to my watch list" link at the top of the page.

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.

About Newsgroups, Newsreaders, and MATLAB Central

What are newsgroups?

The newsgroups are a worldwide forum that is open to everyone. Newsgroups are used to discuss a huge range of topics, make announcements, and trade files.

Discussions are threaded, or grouped in a way that allows you to read a posted message and all of its replies in chronological order. This makes it easy to follow the thread of the conversation, and to see what’s already been said before you post your own reply or make a new posting.

Newsgroup content is distributed by servers hosted by various organizations on the Internet. Messages are exchanged and managed using open-standard protocols. No single entity “owns” the newsgroups.

There are thousands of newsgroups, each addressing a single topic or area of interest. The MATLAB Central Newsreader posts and displays messages in the comp.soft-sys.matlab newsgroup.

How do I read or post to the newsgroups?

MATLAB Central

You can use the integrated newsreader at the MATLAB Central website to read and post messages in this newsgroup. MATLAB Central is hosted by MathWorks.

Messages posted through the MATLAB Central Newsreader are seen by everyone using the newsgroups, regardless of how they access the newsgroups. There are several advantages to using MATLAB Central.

Use the Email Address of Your Choice
The MATLAB Central Newsreader allows you to define an alternative email address as your posting address, avoiding clutter in your primary mailbox and reducing spam.

Spam Control
Most newsgroup spam is filtered out by the MATLAB Central Newsreader.

Tagging
Messages can be tagged with a relevant label by any signed-in user. Tags can be used as keywords to find particular files of interest, or as a way to categorize your bookmarked postings. You may choose to allow others to view your tags, and you can view or search others’ tags as well as those of the community at large. Tagging provides a way to see both the big trends and the smaller, more obscure ideas and applications.

Watch lists
Setting up watch lists allows you to be notified of updates made to postings selected by author, thread, or any search variable. Your watch list notifications can be sent by email (daily digest or immediate), displayed in My Newsreader, or sent via RSS feed.

Other ways to access the newsgroups

Use a newsreader through your school, employer, or internet service provider

Pay for newsgroup access from a commercial provider

Use Google Groups

Mathforum.org provides a newsreader with access to the comp.soft sys.matlab newsgroup