Loren on the Art of MATLAB

Turn ideas into MATLAB

Computational Geometry in MATLAB R2009a, Part II 37

Posted by Loren Shure,

I am pleased to welcome back Damian Sheehy for the sequel to his earlier Computational Geometry post.

 

Contents

What's your problem! ?

Whenever I have a question about MATLAB that I can't figure out from the help or the doc, I only need to stick my head out
of my office and ask my neighbor. Sometimes I wonder how I would get answers to my most challenging how-do-I questions if
I were an isolated MATLAB user tucked away in a lab. I guess the chances are I would get them from the same source, but the
medium would be via the MATLAB newsgroup (comp.soft-sys.matlab).

When I was new to MATLAB I would browse the newsgroup to learn. I found that learning from other people's questions was even
more powerful than learning from their mistakes. As my knowledge of MATLAB has improved, I browse to get an insight into the
usecases and usability problems in my subject area and try to help out if I can. There aren't many posts on computational
geometry, but I have clicked on quite a few griddata subject lines and I was able to learn that for some users, griddata didn't quite "click". Now this wasn't because they were slow; it was because griddata was slow or deficient in some respect.

Despite its shortcomings, griddata is fairly popular, but it does trip up users from time to time. A new user can make a quick pass through the help example
and they are up and running. They use it to interpolate their own data, get good results, and they are hooked. In fact, it
may work so well in black-box mode that users don't concern themselves about the internals; that is, until it starts returning
NaNs, or the interpolated surface looks bad close to the boundary, or maybe it looks bad overall. Then there's the user who fully
understands what's going on inside the black box, but lacks the functionality to do the job efficiently. This user shares
my frustrations; the black box is too confined and needs to be reusable and flexible so it can be applied to general interpolation
scenarios.

griddata Produces "Funny" Results

griddata is useful, but it's not magic. You cannot throw in any old data and expect to get back a nice looking smooth surface like
peaks. Like everything else in MATLAB it doesn't work like that. Sometimes our lack of understanding of what goes on in the inside
distorts our perception of reality and our hopes to get lucky become artificially elevated. It's tempting to go down the path
of diagnosing unexpected results from griddata, but I will have to leave that battle for another day. Besides, if you are reading this blog you are probably a savvy MATLAB
user who is already aware of the usual pitfalls. But, I would like to provide one tip; if 2D scattered data interpolation
produces unexpected results, create a Delaunay triangulation of the (x, y) locations and plot the triangulation. Ideally we would like to see triangles that are close to equilateral. If you see sliver
shaped triangles then you cannot expect to get good results in that neighborhood. This example shows what I mean. I will use
two new features from R2009a; DelaunayTri and TriScatteredInterp, but you can use the functions delaunay and griddata if you are running an earlier version.

I will create a distribution of points (red stars) that would produce concave boundaries if we were to join them like a dot-to-dot
puzzle. I will then create a Delaunay triangulation from these points and plot it (blue triangles).

t= 0.4*pi:0.02:0.6*pi;
x1 = cos(t)';
y1 = sin(t)'-1.02;
x2 = x1;
y2 = y1*(-1);
x3 = linspace(-0.3,0.3,16)';
y3 = zeros(16,1);
x = [x1;x2;x3];
y = [y1;y2;y3];
dt = DelaunayTri(x,y);
plot(x,y, '*r')
axis equal
hold on
triplot(dt)
plot(x1,y1,'-r')
plot(x2,y2,'-r')
hold off

The triangles within the red boundaries are relatively well shaped; they are constructed from points that are in close proximity.
We would expect the interpolation to work well in this region. Outside the red boundary, the triangles are sliver-like and
connect points that are remote from each other. We do not have sufficient sampling in here to accurately capture the surface;
we would expect the results to be poor in these regions.

Let's see how the interpolation would look if we lifted the points to the surface of a paraboloid.

z = x.^2 + y.^2;
F = TriScatteredInterp(x,y,z);
[xi,yi] = meshgrid(-0.3:.02:0.3, -0.0688:0.01:0.0688);
zi = F(xi,yi);
mesh(xi,yi,zi)
xlabel('Interpolated surface', 'fontweight','b');

Now let's see what the actual surface looks like in this domain.

zi = xi.^2 + yi.^2;
mesh(xi,yi,zi)
xlabel('Actual surface', 'fontweight','b');

OK, you get the idea. You can see the parts of the interpolated surface that correspond to the sliver triangles. In 3D, visual
inspection of the triangulation gets a bit trickier, but looking at the point distribution can often help illustrate potential
problems.

Scattered Data Interpolation in the Real-World

Have you ever taken a vacation overseas and experienced the per-transaction exchange-rate penalty for credit card purchases?
Oh yeah, you're on your way out the door of the gift shop with your bag of goodies, you realize you forgot someone and . .
. . catchinggg catchinggg; you're hit with the exchange rate penalty again. But hey, it's a vacation, forget it. Besides,
I'm not Mr Savvy Traveler. Now when you're back at work in front of MATLAB, per-transaction penalties have a completely different
meaning. We have to be savvy to avoid them, but in reality we would simply prefer not to pay them.

The scattered data interpolation tools in MATLAB are very useful for analyzing scientific data. When we collect sample-data
from measurements, we often leverage the same equipment to gather various metrics of interest and sample over a period of
time. Collecting weather data is a good example. We typically collect temperature, rainfall, and snowfall etc, at fixed locations
across the country and we do this continuously. When we import this data into MATLAB we get to do the interesting part – analyze
it. In this scenario we could use interpolation to compute weather data at any (longitude, latitude) location in the country.
Assuming we have access to weather archives we could also make historical comparisons. For example, we could compute the total
rainfall at (longitude, latitude) for years 2008, 2007, 2006 etc, and likewise we could do the same for total snowfall. If
we were to use griddata to perform the interpolation, then catchinggg catchinggg, we would pay the unnecessary overhead of computing the underlying
Delaunay triangulation each time a query is evaluated.

For example, suppose we have the following annual rainfall data for 2006-2008:

     (xlon, ylat, annualRainfall06)
     (xlon, ylat, annualRainfall07)
     (xlon, ylat, annualRainfall08)

xlon, ylat, and annualRainfall0* are vectors that represent the sample locations and the corresponding annual rainfall at these locations. We can use griddata to query by interpolation the annual rainfall zq at location (xq, yq).

     zq06 = griddata(xlon, ylat, annualRainfall06, xq, yq)
     zq07 = griddata(xlon, ylat, annualRainfall07, xq, yq)
     zq08 = griddata(xlon, ylat, annualRainfall08, xq, yq)

Note that a Delaunay triangulation of the points (xlon, ylat) is computed each time the interpolation is evaluated at (xq, yq). The same triangulation could have been reused for the 07 and 08 queries, since we only changed the values at each location.
In this little illustrative example the penalty may not be significant because the number of samples is likely to be in the
hundreds or less. However, when we interpolate large data sets this represents a serious transaction penalty. TriScatteredInterp, the new scattered data interpolation tool released in R2009a, allows us to avoid this overhead. We can use it to construct
the underlying triangulation once and make repeated low-cost queries. In addition, for fixed sample locations the sample values
and the interpolation method can be changed independently without re-triangulation. Repeating the previous illustration using
TriScatteredInterp we would get the following:

     F = TriScatteredInterp(xlon, ylat, annualRainfall06)
     zq06 = F(xq, yq)
     F.V = annualRainfall07;
     zq07 = F(xq, yq)
     F.V = annualRainfall08;
     zq08 = F(xq, yq)

We can also change the interpolation method on-the-fly; this switches to the new Natural Neighbor interpolation method.

     F.Method = 'natural';

Again, no transaction penalty for switching the method. To evaluate using natural neighbor interpolation we do the following:

     zq08 = F(xq, yq)

As you can see, the interpolant works like a function handle making the calling syntax more intuitive.

What's the Bottom Line?

The new scattered data interpolation function, TriScatteredInterp, allows you to interpolate large data sets very efficiently. Because the triangulation is cached within the black-box (the
interpolant) you do not need to pay the overhead of recomputing the triangulation each time you evaluate, change the interpolation
method, or change the values at the sample locations.

What's more, as in the case of DelaunayTri, this interpolation tool uses the CGAL Delaunay triangulations which are robust against numerical problems. The triangulation algorithms are also faster and more
memory efficient. So if you use MATLAB to interpolate large scattered data sets, this is great news for you.

The other bonus feature is the availability of the Natural Neighbor interpolation method – not be confused with Nearest Neighbor
interpolation. Natural neighbor interpolation is a triangulation-based method that has an area-of-influence weighting associated
with each sample point. It is widely used in the area of geostatistics because it is C1 continuous (except at the sample locations)
and it performs well in both clustered and sparse data locations.

For more information on this feature and the new computational geometry features shipped in R2009a check out the overview video, and Delaunay triangulation demo, or follow these links to the documentation

Feedback

If you use MATLAB to perform scattered data interpolation, let me know what you think here. It's the feedback from users like you that prompts us to provide enhancements that are relevant to you.

Get the MATLAB code

Published with MATLAB® 7.8

Note

Comments are closed.

37 CommentsOldest to Newest

Andrew Scott replied on : 2 of 37
In reconstruction of magnetic resonance images (MRI) data is acquired in the spatial frequency domain (k-space). In some cases the data is acquired along spiral or radial paths in k-space and gridding is used to interpolate the data onto a Cartesian grid so a fast Fourier transform can be used to obtain an image. Commonly the gridding procedure involves convolving the data with a gridding kernel (often a Kaiser-Bessel function), compensating for the density of sample points in k-space and compensating for the effect of the convolution on the final image after the FFT. I have, however, managed to get images using the MATLAB griddata functions but don't really understand how the Delaunay triangulation works so wrote my own gridding routines. Do you have a feel as to how these two approaches vary? Might a future MATLAB version incorporate gridding via convolution?
Damian Sheehy replied on : 3 of 37
Hi Andrew, Whenever I come across MRI images in graphics journals and magazines like New Scientist (www.newscientist.com) the vivid colors capture my interest and leave me feeling awestruck by this technology. I often wish I had some background in medical imaging to understand the algorithms, but I don’t so I cannot provide an insight into how these gridding algorithms compare. When I was writing this blog it struck me that very few textbooks cover Delaunay-based interpolation. Journal publications on the subject assume you know, and in case you don’t they usually cite David Watson’s out-of-print book as an introductory reference. Last time I checked on Amazon, a used copy of Watson’s little book - “Contouring: A guide to the analysis and display of spatial data”, Pergamon, 1994 – was selling for over a hundred and forty US dollars. Delaunay-based interpolation is very straightforward. Suppose you have a set of sample points as follows:
t = 0:pi/4:2*pi;
a = t(1:length(t));
x = sin(a)';
y = cos(a)';
x(end) = 0;
y(end) = 0;
plot(x,y,'ok')
axis equal
Which we can triangulate like this:
hold on
tri = delaunay(x,y);
triplot(tri,x,y)
Now, to interpolate a value at some location such as this. . .
plot(0.25,-0.6,'*r')
We simply search for the triangle enclosing the point. To perform linear interpolation we weight the values of each vertex of the enclosing triangle using an area-based weighting. There are various ways to connect the dots to build the triangulation, but the Delaunay algorithm is the best because it favors connecting points that are in close proximity to each other. For example, it will choose the configuration in Figure 1 over the one on Figure 2.

close(gcf)
x = [0 0.5 1 0.5]';
y = [0 -2 0 2]';
tri = [1 2 4; 2 3 4];
triplot(tri,x,y)
axis equal

tri = [1 3 4; 1 2 3];
figure
triplot(tri,x,y)
axis equal

Also, when performing the actual interpolation there are various ways to choose the vertices of influence and the corresponding weights for each vertex, but you get the gist. In 3D the Delaunay triangulation constructs tetrahedra as opposed to triangles, but the theory extends naturally. I hope this gives you some insight into how the Delaunay-based interpolation works. Damian
Damian Sheehy replied on : 4 of 37
Second-last paragraph of my previous entry should have said For example, it will choose the configuration in Figure 2 over the one on Figure 1. Apologies for the confusion.
Aslak Grinsted replied on : 5 of 37
Thanks. this is very useful. Trianglulations are not always the best method for gridding a data set like precipitation. I really wish that matlab would offer more gridding methods: Kriging, radial basis functions, ... Basically i would like to be able to do what Surfer does, only inside matlab. I like gridfit.m on fileexchange, but I use my own surfergriddata.m (also on fileexchange) which is just an activex-wrapper for surfer. I suspect both of those will perform much better on the sliver example.
Damian Sheehy replied on : 6 of 37
Hi Aslak, It’s good to hear from you again! Thank you for your feedback, and I am happy to hear you will find the new interpolation tool useful. Yes, triangulation-based interpolation is one of several useful techniques, and it does have certain limitations. I chose the weather data use-case as a readily-understood example to illustrate the flexibility of the new interpolation tool, rather than advocating the use of triangulation techniques in this application area. But, I’m sure you’ll forgive me for this. More importantly, your enhancement request to support interpolation based on Kriging, and radial basis functions is something we will consider. We have had several requests for these features, and in fact it was requests from users like you that prompted us to support the Natural Neighbor technique we provided in TriScatteredInterp. Damian
Jason Stockmann replied on : 7 of 37
Is there a version of griddata that allows the user to choose from more types of interpolation, or better yet, specify their own interpolation kernel as an input? This would allow the user to use to convolve irregularly-spaced data with a kernel and interpolate onto desired points on a regular grid. This would be very useful for experimenting with different types of kernels for regridding MRI k-space data from spiral, radial, rosette, random, and sundry other trajectories.
Damian Sheehy replied on : 8 of 37
Hi Jason, The interpolation methods provided by GRIDDATA and TriScatteredInterp are triangulation-based. The 'v4' method in GRIDDATA is an exception, and is not carried over to TriScatteredInterp. I am not aware of variants on the file exchange that extend GRIDDATA to address your usecase. Andrew Scott also cited this usecase in the second comment. I have a colleague here at The MathWorks who has a strong background in Medical Imaging and I will discuss this enhancement request with him to gain some insight into how we can better support this workflow. Thank you for your feedback and suggestion! Damian
Adnan replied on : 9 of 37
I´m working with complex data interpolation. The function ´griddata´was able to take the complex data, but 'TriScatteredInterp' doesn´t. For my application, I can not use griddata (due to certain precision issues with qhullmx sub-routine). At the same time,'TriScatteredInterp' won´t take my complex data. Is there a way around? Or can you recommend me some other interpolation function in matlab that would accept complex data?
Damian Sheehy replied on : 10 of 37
Hi Adnan, You just need to interpolate the real and imaginary parts separately and combine the results afterwards. Here's an example: % Create some sample data X = -5 + 10*rand(10,1); Y = -5 + 10*rand(10,1); complexV = complex(rand(10,1), rand(10,1)); % Suppose we want to evaluate the interpolant at query points Qx, and Qy % We will just query at 5 random locations Qx = rand(5,1); Qy = rand(5,1); % Use REAL and IMAG to extract the real and imaginary components. realV = real(complexV); imV = imag(complexV); % Create the interpolant that will interpolate over the real data F = TriScatteredInterp(X,Y,realV); % Interpolate the real data at the query locations Qvreal = F(Qx,Qy); % Now modify the interpolant to interpolate over the imaginary data F.V = imV; % Interpolate the imaginary data at the query locations Qvim = F(Qx,Qy); % The complex result is Qv = complex(Qvreal, Qvim) Best regards, Damian
James replied on : 11 of 37
I understand that Natural Neighbour interpolation may meet problem when the interpolation points are close to the boundary of original data points. Howvever, I tried different original and interpolation datasets (irregular or regular), also used INTERP2 with "v4" and "cubic". For all the datasets, Natural Neighbour interpolation with TriScatteredInterp showed similar accuracy to linear interpolation and much worse than INTERP2 with "v4" and "cubic". Could you please provide me an example matlab code which compares Natural Neighbour interpolation with "spline" or "cubic" and shows Natural Neighbour interpolation is superior?
Joris Timmermans replied on : 12 of 37
I am working with Remote sensing data that needs to be resampled to a rectangular lat lon grid. As the image has a non-normal projection, the usual imtransform options do not work for me. I therefore tried griddata (and more recently the TriScatteredInterp). Why is the "TriScatteredInterp" method slower for large datasets than "griddata"?
Damian Sheehy replied on : 13 of 37
James, Each interpolation method has its own relative merits. There are various factors to consider; performance, accuracy of fit, continuity, distribution of sample data, etc. MATLAB provides the methods and the user decides which method best meets their needs in the context of the application at hand. I suggest you review some recent literature that deals with a comparison of the various methods in the area of scattered data interpolation. Damian
Damian Sheehy replied on : 14 of 37
Joris, If you are using release R2009a, TriScatteredInterp may exhibit slower performance if you have a large sample data set that lies on a regular grid or parallel scan lines. We made improvements to the performance in release R2009b and R2010a that addresses this issue. In the R2010a release, TriScatteredInterp is generally 3 to 4 times faster than the GRIDDATA function of R2009a or before, and should not be appreciably slower under any circumstances. The R2010a prerelease is currently available. Damian
Joris Timmermans replied on : 15 of 37
Thank you very much for you answer. I have asked our IT departement to install the new version on the netwerk.
Abhishek replied on : 16 of 37
I am working with a problem that requires me to construct the convex hull of a multidimensional scattered data-set and then compute the intersection of that with a (hyper)plane of type x1 = constant in n-dim, where I know the numerical value of this constant. I realize that the intersection is a (n-1) dim convex polytope (for example, you get a convex polygon when you cut a 3d hull by x = constant plane) but am wondering how to compute it. Notice that, this convex region-of-intersection is different than the projection of the hull onto the specified plane.
Damian Sheehy replied on : 17 of 37
Hi Abhishek, You did not mention how many dimensions you are working in, so let me start by saying the convex hull algorithm is designed to handle small to medium dimensions and by that I mean less than 10. For a large point set you may encounter the Curse of Dimensionality in dimensions lower than 10. Once you have computed the convex hull extract the boundary edges. For each edge test the end points against the plane. If the end points lie on opposite sides of the plane, then the edge intersects the plane; compute the intersection point and add it to the set X. Finally, compute the convex hull of the 2D projected point set Xproj. Consider the 3D case, if your plane is defined by a point (xi, yi, zi) and normal (a, b, c) substitute these values into the equation ax + by + cz + d = 0 to compute d. To test a query point, (xq, yq, zq), against the plane substitute (xq, yq, zq) into the equation; if the result is positive the point lies on the positive side of the plane, and if it’s negative it lies on the opposite side. To compute the intersection point, refer to the following; http://en.wikipedia.org/wiki/Line-plane_intersection Everything should extend naturally to higher dimensions. I hope this helps to get you on the right track. Damian
Karen Bemis replied on : 18 of 37
will TriScatteredInterp do true 3D interpolation, that is, will it interpolate a set of scattered points value@(x,y,x) to a uniform mesh? I'm starting out with values at (r,theta,phi) or something sort of similar. I calculate their x,y,z positions but then want to interpolate to a uniform spacing in x,y,z for visualization and analysis. Maybe I'm missing something but I can't figure out how to do this with TriScatteredInterp. Actually, I'm not sure griddata could do it (we've been using a in house nearest neighbor interpolation routine).
Damian Sheehy replied on : 19 of 37
Hi Karen, Yes, TriScatteredInterp will interpolate scattered data in 3D space and allow you to evaluate the interpolant at gridpoint locations. In fact, you can evaluate the interpolant at any (x,y,z) location within the convex hull of your sample data points. Perhaps you are trying to perform interpolation on the surface of a sphere? Are you trying to convert scattered data on the surface of a sphere into a grid defined in Longitude-Latitude space? If so this is a different problem because it amounts to interpolating a 2D manifold embedded in 3D space. This is type of problem is not supported by griddata or TriScatteredInterp. Google "Interpolation on the surface of a sphere" for more information on this particular subject. Damian
Joris Timmermans replied on : 20 of 37
If I use TriScatteredInterp for resampling imagery, I bump into some problems. -------------------------------------------------------- Using Lon(518400x1),(min(Lon(:))~-5.14, max(Lon(:))=08.99 Lat(518400x1),(min(Lat(:))~40.00, max(Lat(:))=60.00 Values_t0(518400x1) F= TriScatteredInterp(Lon,Lat,Values_t0,'nearest'); imagesc(F(Loni,Lati)); -------------------------------------------------------- *First it gives my a warning: "Warning: Duplicate data points have been detected and removed - corresponding values have been averaged." *But when I plot the data it provides me with the image that I want. However I want to resample multiple timeframes using the same triangulation (reducing the per-transaction exchange-rate penalty): -------------------------------------------------------- F.V = Values_t1(:); imagesc(F(Loni,Lati)); -------------------------------------------------------- *It gives me an error: "??? Error using ==> subsref The interpolant is in an invalid state.The number of data point locations should equal the number of data point values." *After inspection it appears that the TriScatteredInterp provides me with size(F.V)= 516952x1 instead of 518400x1. How can I circumvent this problem with kind regards Joris Timmermans
Damian Sheehy replied on : 21 of 37
Joris, Let's follow a simple example that contains duplicates and see what's going on. . .
% 1 - Create sample point locations
X = rand(20,1);
Y = rand(20,1);

% 2 - Duplicate some point locations
X(14) = X(4);
Y(14) = Y(4);
X(16) = X(6);
Y(16) = Y(6);

% 3 - Create some corresponding values
V = X.^2 + Y.^2;

% 4 - Now create the interpolant
F = TriScatteredInterp(X,Y,V);

% We get the following warning:
Warning: Duplicate data points have been detected and removed - 
corresponding values have been averaged.
When we check the Interpolant we get the following:
 F

F = 

  TriScatteredInterp

  Properties:
         X: [18x2 double]
         V: [18x1 double]
    Method: 'linear'

  Methods
We can evaluate the interpolant and it works as expected
Vq = F(0.1, 0.2)

Vq =

    0.0859
At this point the interpolant is in a valid state. It does not "remember" that it had to remove points 14 and 16 to get to this state - why should it? Now we wish to change the values - we will attempt to assign 20 values to an interpolant that has 18 data points.
V = 3*(X.^2) + Y.^2;
F.V = V
Now check the interpolant
F

F = 

  TriScatteredInterp

  Properties:
         X: [18x2 double]
         V: [20x1 double]
    Method: 'linear'
Now try to evaluate the interpolant
Vq = F(0.1, 0.2);

??? Error 
The interpolant is in an invalid state.
The number of data point locations should equal the number of data point values.
In the first scenario MATLAB warned me I had duplicate locations, but it did the right thing and I ignored the warning. In the second scenario the unheeded warning gave rise to problems because MATLAB can't read my mind and doesn't understand my intent. The correct approach is to heed the warning and remove duplicates first, and then work with a unique data from there on. So back up a step and follow that approach. We can use the unique function to eliminate duplicates.

V = X.^2 + Y.^2;
[~, I, ~] = unique([X Y],'first','rows');

I = sort(I);
X = X(I);
Y = Y(I);
V = V(I);
F = TriScatteredInterp(X,Y,V) % No warnings
Now assign new values.
V = 3*(X.^2) + Y.^2;
F.V = V
Now evaluate the interpolant and everything is OK
Vq = F(0.1, 0.2); % OK
Ben replied on : 22 of 37
Nice post. I had feared that any interpolation scheme might do poorly in a large concave region, so it was nice to see this confirmed. I have a set of 2D experimental data where most of the interior of the domain is filled in (a few holes here and there) but the edges are especially jagged. I would rather not interpolate at all in the concave regions along the edges, or in the interior holes, if the hole is big enough. Is there an interpolation routine where I can force it to only interpolate from points within a specified radius? If it does not have enough points within the radius, then I would like the routine to return NaN. This way I would hopefully avoid the slivers in the DelaunayTri example above. At the moment all I can think of is to somehow find the outside edges of the 2D domain, do a constrained DelaunayTriangulation, and query the inOutStatus to locate the triangles within the domain. Is there a better way? It's funny. It seems like everyone is interested in getting a better extrapolation/interpolation and there are lots of options for that. But, I'm having a lot of trouble finding tools for the more conservative approach of only interpolating if you have points close enough by.
Ben replied on : 23 of 37
Ok, I just found all the other TriRep functions. Some combination of circumcenters and neighbors should do the trick for me. You can disregard my last comment.
Huseyin replied on : 24 of 37
The option to avoid multiple triangulation seems useful. However I am getting some strange results.
size(Xp)
tic
F = TriScatteredInterp(Xp,Yp,Hsig);
toc
Hsig2=Hsig;
tic
F.V=Hsig2;
toc
OUTPUT: ans = 47815 1 Elapsed time is 0.408701 seconds. Elapsed time is 0.471277 seconds. Updating V takes more time than TriScatteredInterp. Can someone explain this?
Mike replied on : 25 of 37
I'm having similar results as Huseyin -- updating F.V takes slightly longer than creating F with TriScatteredInterp. It's still faster than repeated calls to griddata, but it sounds like it should be quicker, based on the description above about "repeated low-cost queries".
Damian Sheehy replied on : 26 of 37
Huseyin, Mike, The optimization related to the editing of the properties is only realized when the code resides in a file. This is the general mode of execution for development applications. Drop the following function in a file named ReplaceVTest.m
function ReplaceVTest()
  x = rand(1000000,1)*4-2;
  y = rand(1000000,1)*4-2;
  z = x.*exp(-x.^2-y.^2);
  tic; F = TriScatteredInterp(x,y,z); toc
  tic; F.V = 2*z; toc
When you execute the test you get the following > ReplaceVTest Elapsed time is 16.406036 seconds. % Time to create the interpolant Elapsed time is 0.014721 seconds. % Time to replace the values
Huseyin replied on : 27 of 37
Thanks Damien. I don't get it but it works. And I guess you mean 'when the code resides in a "function".'
Damian Sheehy replied on : 28 of 37
Huseyin, That’s correct, the code must reside in a function to achieve optimal performance. When MATLAB executes a program that is composed of functions that reside in files, it has a complete picture of the execution of the code; this allows MATLAB to optimize for performance. When you type the code at the prompt MATLAB cannot anticipate what you are going to type next so it cannot perform the same level of optimization. Developing applications through the creation of reusable functions is general and recommended practice, and MATLAB will optimize the performance in this setting.
Huseyin replied on : 29 of 37
Damian, this is interesting. Are there other functions that benefit from function optimization as much as this? How do we predict if a function will be faster inside a function rather than a script?
Mike replied on : 30 of 37
My code was in a .m file, but it was not a function. I added only a "function name()" line to the top, and it cut execution time from 7 hrs to 7 mins. Thanks, Damian. Is such a performance improvement attributable to functions being fully parsed before execution, while scripts are not?
Damian Sheehy replied on : 31 of 37
Huseyin, Mike, It is difficult to quantify the performance improvements that can be gained by modularizing your code into functions. It depends on several factors including the functions you are calling, the size of the data, and level of iteration involved, etc. The large performance gains that Mike observed are probably unique to the computational geometry classes in the context of problem size and level of iteration. MATLAB code is designed to execute with optimal performance when structured into functions, so you do not need to be concerned about the internals if you adopt this recommended approach when developing MATLAB applications.
Mike replied on : 32 of 37
Hi, I'm trying to use TriScatteredInterp, but there's a catch. The x,y points apparently has some repetition. When I run TriScatteredInterp, I get the following message: "Warning: Duplicate data points have been detected and removed - corresponding values have been averaged." That's all fine and dandy. (technically though, the message shows up even if there are no values to average - i.e. if I just introduced x,y data by F.X = XY) But the problem is that I cannot reuse this interpolation F to run over the same x,y points with a different set of z-data, because the other z-data defined over the same x,y will not have the redundant points removed. And that only works if the z-data and x,y are introduced together via the function - F = TriScatteredInterp(X, V). First getting an empty F, then setting F.X = XY, then F.V = data will not work. If I attempt to use F with other z-data, I get: "??? Error using ==> subsref The interpolant is in an invalid state. The number of data point locations should equal the number of data point values." In this case, it looks like I have to create a new interpolant for every set of z-data, which seems highly inefficient, and is only due to a trivial problem of not knowing how to remove the redundant points in the same way. So I was wondering if you had a solution for this. I could try to manually remove redundant points beforehand, but then I have the added overhead of remembering some indexing matrix for removing redundant points and applying it to any future z-data that may be used, and also averaging the redundant values. I would suggest that it would be a lot cleaner to internalize this somehow. Since the interpolant is designed to be re-used for different z-data over the same x,y points, it should remember if it had to pre-process the x,y points. Then, any z-data introduced by F.V = data should undergo the same pre-processing, so that the user doesn't have to worry about nitty gritty details like this. Maybe I'm just using this wrong and there is an easier way, in which case please let me know!
Damian Sheehy replied on : 33 of 37
Mike, In item #21 above, I outlined the approach to handling duplicate data points in the scenario where the interpolant is to be reused with different values V. This suggests removing the duplicates prior to performing the interpolation and the steps to do this are provided. TriScatteredInterp allows X and V to behave independently of each other. This flexibility allows the underlying triangulation to be reused when V is changed and this is very desirable because a retriangulation of X incurs an overhead. In your duplicate point scenario X and V are not strictly independent of each other and the design cannot handle both situations faithfully. TriScatteredInterp cannot and should not remember which points were merged when the interpolant was created, because when the interpolant is subsequently edited we cannot assume intent. The user may wish to replace V and retain the same locations X, OR the user may wish to replace V and then replace X with a completely different set of locations. A more conservative design would restrict you to editing both X and V simultaneously; e.g. F.SetData(newX,newV), and this has the same overhead as GRIDDATA. Thanks for taking the time to provide feedback, we will use it to improve the documentation for this feature.
Darin Williams replied on : 34 of 37
Your blog makes some great points about the perils of triangulated interpolation. Suppose that "suspect" boundary triangles are already identified, and one would like to treat them as though they were outside the hull for interpolation, returning points withing them as nans. How could this be done with triscatteredInterp? Note thaat this, effectively imposes a non-convex boundary on the interpolation.
Damian Sheehy replied on : 35 of 37
Darin, This is a good question on a usecase I can see has practical relevance. When we design a feature like this we try to strike a balance between functional ease-of-use and complexity associated with the more advanced bells and whistles. Your usecase falls into the advanced category and while TriScatteredInterp doesn't readily provide a "canned" solution you can do this with a couple of extra steps. We will reuse the above example that illustrated boundary slivers and create a DelaunayTri from the data as before.
t= 0.4*pi:0.02:0.6*pi;
x1 = cos(t)';
y1 = sin(t)'-1.02;
x2 = x1;
y2 = y1*(-1);
x3 = linspace(-0.3,0.3,16)';
y3 = zeros(16,1);
x = [x1;x2;x3];
y = [y1;y2;y3];
dt = DelaunayTri(x,y);
We can use the aspect ratio (AR) of the triangle as a rejection criterion; the AR being defined as the ratio of the incircle radius to circumcircle radius. I did some experimenting and found 0.01 was a good threshold for this example.
[~, Ric] = dt.incenters();
[~, Rcc] = dt.circumcenters();
AR = Ric./Rcc;
Now find the indices of the admissible triangles and visually inspect.
Tadmissible = find(AR > 0.01);
trisurf(dt(Tadmissible  ,:),x,y,zeros(size(x)))
Assign NaNs to the points that are not within our region of interest. We will use the same query data as before and use the pointLocation method provided by DelaunayTri to find the triangles enclosing the points.
[xi,yi] = meshgrid(-0.3:.02:0.3, -0.0688:0.01:0.0688);
gridsize = size(xi);
xi = xi(:);
yi = yi(:);
% Initialize the corresponding values using NaNs and later replace the valid entries.
zi = NaN(size(xi,1),1);
Tenclosing = dt.pointLocation(xi,yi);
validxiyi = ismember(Tenclosing,Tadmissible );
You can reuse the DelaunayTri by passing it to TriScatteredInterp - don't waste anything you can reuse!
z = x.^2 + y.^2;
F = TriScatteredInterp(dt,z);
%
% Evaluate the interpolant at the valid query locations
zi(validxiyi) = F(xi(validxiyi), yi(validxiyi))
% Reshape the data so we can call mesh
xi = reshape(xi,gridsize );
yi = reshape(yi,gridsize );
zi = reshape(zi,gridsize );
mesh(xi,yi,zi)
xlabel('Interpolated surface', 'fontweight','b');
If linear interpolation suffices, you could use DelaunayTri to perform the interpolation. The pointLocation method returns the Barycentric coordinates bc which are used to weight the values at each vertex of the enclosing triangle.
dt = DelaunayTri(x,y);
z = x.^2 + y.^2;
[ti bc] = dt.pointLocation(xi,yi);
T = dt.Triangulation;
zi = z(T(ti,1)).*bc(:,1)+z(T(ti,2)).*bc(:,2)+z(T(ti,3)).*bc(:,3)
The user documentation for Scattered Data Interpolation has been revised recently and may help clarify the latter approach.
Darin replied on : 36 of 37
Poblem associating X/Y/Z poitns after triangulating in XY via DelaunayTri. If the X/Y contain duplicate points they are removed. You no longer know which Z goes with which XY, or what DT.X goes with the input XY. Is DelaunayTri's criteria for removal EXACT equality of X and Y? Is there an easier way than pre-screening to remove matches? (Seems to defeat some of the improvment over delaunay.) How about adding optional output parameters to descramble this?
[DT,ind] = DelaunayTri(X,Y);  %??
where ind is a list of the DT.X output index receiving each XY input point (for use in accumarray)?
Damian Sheehy replied on : 37 of 37
Darin, When you call DELAUNAY with a set of points that contain duplicates you will discover some data points are not referenced by the triangulation. If you were to use that triangulation to develop an algorithm, this could cause problems. Adding a method to DelaunayTri to provide information about duplicates may seem like a good idea at first. But in terms of interface design we are really concerned about the separation of responsibility. DelaunayTri is for triangulating data and querying that triangulation. Eliminating duplicates is the responsibility of another function. Besides, you would still need to intervene to get what you want. For example, when eliminating duplicates you may wish to average the corresponding Z values.