Comments and Ratings

Author comment: v24 was posted to resolve a licence conflict. The Mathworks no longer allows code to be released under a GNU GPL, so this has been removed for v24.

There are some (minor) revisions to the code, based on user feedback. Specifically, v24 may be slightly slower than v23, but should generally produce higher quality meshes, especially for complex geometries.

Comments/feedback is always welcome - d_engwirda@hotmail.com

I would like to have more information of the theory behind this. What it means that it is "quadratic"? Should it be exact for 2nd. degree functions? I tried this
X=[1:10 1:10];Y=[zeros(1,10) ones(1,10)];%Data points
f=100*(X').^2; %data points
tri=delaunay(X,Y);%Triangulation
triplot(tri,X,Y) %to interpret the case
x=linspace(0,11,300);y=repmat(0.5,1,300);
fC=tinterp([X' Y'],tri,f,x,y,'quadratic');
plot(diff(fC));%should be linear if the interpolant is exact for 2nd degree

Apparently from this example the interpolant is C0 (i.e. continuous but not differentiable), but it is a little bit misleading what is meant by "quadratic".

I tried to contact the author before posting, but apparently changed addres. Thank you anyway for his contribution.

There's actually a very important point here. Because of how Matlab stores a sparse matrix (see Tim Davis' book, or read the help on, say, mxGetIc), and because a major bottleneck for these computations is loading data into the CPU's cache, it is MUCH faster to store a matrix by row, instead of column, when doing matrix-vector multiples.

In Matlab, this is easy to remedy: simply take the transpose.

So, this cool trick can save you time. If you want to do this repeatedly:
>> y = A*b
instead, do this:
>> At = A.'; % a one-time cost
>> y = At.'b

But, this trick only works on newer versions of Matlab. On, say, R2006b, a call like
>> y = At.'b
will actually calculate the transpose of At, so this is very slow!

On older versions of Matlab, it's likely that this SMVP code will really help, but it might not be as helpful on R2008, for example.

On linux, non-multithreaded R2006a, I found that SMVP took about 60% the time of Matlab's own multiply. So, good work!

Thank You,
 it was very helpfull, impressive code.

If I can suggest an improvement:
a mex version will be the top of performance.

It is a good job, but I still find some problems when detecting a point of a polygon lies on the polygon edge. My test run in Matlab R2006a. Suppose node be the points of a polygon,
"[in,on] = inpoly(node,node);" can get right results, but after running
"[in,on] = inpoly(node(1,:),node);", on is false.

Very fast, very accurate. Good jobb, Darren!

inpoly provides a very large speed increase for large polygon problems!

Unfortunately, I too noticed a bug.

I am checking a very large rectangle grid's points to see if they are in a "massively concave polygon" -- think the outline of a "robot with arms, legs, etc."

The inpoly algorithm "incorrectly" *ADDS* a "shock of hair" to the "robot" -- obviously I am speaking metaphorically here -- and hence I am stuck using the *MUCH SLOWER* "inpolygon".

Has anyone reported an error like this to you?

Is there anything I can do?

Is there some middle ground between the two -- for example, by and large, I am dealing with concave polygons that just have an outside boundary -- no hole or anything.

Small bug (as noted below) fixed.

Users don't beware, inpoly should work in all cases.

Further bugs can be emailed if necessary...

It's an incredible speed up compared with the poor inpolygon function delivered with MATLAB! Must be O(M*log(N)) ;-)

Works but offers almost very minor 8% speedup as of Matlab 2006.
I multiplied 700,000 x 700,000 matrix with 2,000,000 nonzero entries times a 700,000 vector and got only 8% faster than Matlab 2006. Not worth using unless there is more speedup in other settings.

2) The routine fails to check that the matrix allociation succeeded (may run out of memory). Need to add a NULL check to

plhs[0] = mxCreateDoubleMatrix(mxGetM(prhs[0]), 1, mxREAL);

3) Change if (nlhs!=1) to if (nlhs > 1)

Summary: works but very limited speedup.

An alternative to griddata - similar in speed when you add in the time for the triangulation.

The advantage of this code is if you already have a triangulation, then this saves the time to regenerate it. It also allows you to use a custom, non-delaunay triangulation - not an option at all for griddata.

A disadvantage of this code is it forces you to do the triangulation to call it. Its an unnecessary step for many people, so I'd probably be tempted to allow the user to not supply the triangulation, so if t is empty it could just call delaunay.

This code offers properties similar to griddata. It will not extrapolate beyond the convex hull of the data, Extrapolation is risky business of course, so thats not necessarily bad.

Good help, examples. The quadratic interpolant is an interesting idea.

It's really faster. Unfortunately, it does
not support complex matrices.

I had no problems using this file and it is around 5 times quicker than matlab's A*b for the problems I tried.

Very good contribution. I wish the author could extend the tricontour function to handle 3D triangle meshes as well.

QUOTE
The reason that cnect is defined on its own is so that the domain can be multiply connected (polygon with "islands").
QUOTE
please, note that i said at run-time, i still feel that you should make the third arg optional
us

While I prefer the edge list implementation this code uses as opposed to Matlab's polygon, it would be easy enough to generate the connectivity assuming consecutive points on the polygon as us points out.
Regardless, this code is indeed fast and nice.