Surface Fitting using gridfit

Model 2-d surfaces from scattered data

Description

Those wishing to model a surface from data in the
form of z(x,y) from scattered or semi-scattered
data have had few options in matlab - mainly
griddata.

Griddata is a valuable tool for interpolation of
scattered data. However it fails when there are
replicates or when the data has many collinear
points. Griddata is also unable to extrapolate
beyond the convex hull of the data unless the 'v4'
option is used, which is slow.

Gridfit solves all of these problems, although it
is not an interpolant. It builds a surface over a
complete lattice, extrapolating smoothly into the
corners. You have control of the amount of
smoothing done, as well as interpolation methods,
which solver to use, etc.

This release allows the user to solve much larger problems using a new tiling option. There is essentially no limit on the size of the suface one builds now, as long as you have dense enough data and enough memory to store the final gridded surface.

Example uses are found in the file gridfit_demo.m,
as well as comparisons to griddata on the same
surfaces.

• MATLAB release: MATLAB 7.0.1 (R14SP1)
• Other requirements: Gridfit should be accessible to older versions of matlab, as long as symamd and colamd are available. Older releases could use one of the iterative solvers, such as symmlq.

[Michael Robbins]

Neat picture!

[Urs (us) Schwarz]

simply said: just another extremely useful member of the d'erriconian-draconian family of must-haves if you're professionally serious about nd-data investigation / reconstruction / interpolation / visualization, which smoothly unites with its now-famous siblings (consolidator, inpaint-nans) into a very nice toolbox for those everyday cheating-with-data endeavors...
the code is very clean and the profiler reveals no unnecessary hot-spot...
HOWEVER, in contrast, the help bit is just unbelievably clumsy and almost indigestible without sever anxiolytic medication - JOHN, JOHN! this part requires some serious overhaul and sprucin'-up...
altogether, very well done - and thank you for this snippet
us

[Chaowen Liang]

This is exactly what I am looking for... Thanks you!!! Save me a lot of time to write one...

[Alberto Baraldi]

I was struggling against griddata since days, beacause of its problems when extrapolating...this is absolutely great!
Thank you for your great help to the Matlab community. We love you

[Kong Zuor]

Well done John, I have been looking for this while ago and I think it is a must have code. Thank you very much for being so helpful to mathlab community

[CASTRY Alen]

Thank you...

[Julio Oliveira]

Great ! Thank you!

[David Barker]

This is an excellent program. It makes it easy to fit surfaces to 3D data.

[Kevin Denis]

 

[T O]

Excellent Work!

[Wen Soong]

Does a great job!

[Aslak Grinsted]

Should be included into the default distribution of matlab! perhaps somehow integrated into griddata?

[thank you]

Marvelous!

[Andy Lulham]

 

[Håvard Torpe]

Works like a charm right from the moment you start using it.

However: I'd like gridfit() NOT to extrapolate values for gridnodes that has no datapoint in the vicinity. Any ideas?

[John D'Errico]

I'll invite Håvard Torpe to contact me via e-mail to discuss this idea. But while I COULD allow the user to turn off extrapolation in gridfit, this will likely cause singularities - serious numerical problems, for most users who might do so.

[Christopher Schwalm]

Excellent program and documentation.

A question: the help states "Gridfit is not an interpolant." Why is that? It seems to be perfectly able to be used as such, given the outputs and how they were generated. One use I see in my work is to use this function to interpolate where repeated values are present. I've yet to find another Matlab solution that works as well as this does.

[Wolfgang Schwanghart]

Excellent! This is exactly what I was searching for. Thank you very much for this great contribution.

[Phani Ivatury]

Excellent!!!

[Thomas Clark]

Excellent function, John - this has been extremely useful to me.

The ability to turn off extrapolation beyond the bounds of the data would be very useful, though.

A note for the inexperienced user:
If you have repetitions (multiple points with the same (x,y) but different (z)), this algorithm will handle where griddata might fall. In such a case, think hard about why you have these repetitions; it may be that volumetric plotting is more appropriate than surfaces. If so, you could get a good-looking surface out of this algorithm; which isn't totally meaningful.

That's not a fault with this algorithm; it just means that you have to think hard about why you are getting repeated points in Z: If it's acceptable to average these out to a surface, this algorithm is for you!

[Carlos Adrián Vargas Aguilera]

Very very good John! Take good care of Amy!

[E Farhi]

This is simply great, with all we can expect for a replacement to griddata (which fails too often).

[Zhijun Wang]

This is very very good work!

[Ken Lin]

Thanks for your programs. It solves my problem now:)
by the way, would you like to work out some way to interpolate the data just like the surfer does

[David Sprinkle]

Great tool and insightful docs! Gridfit succeeds where griddata fails; that is, with noisy & ill-spaced data such as exist in the real world. Thanks!

[gabriele flauti]

 

[Juan Perez]

 

[Hendawi Mohamed]

I just write to say..."Thank you"!

PS. but I add: It would be perfect if there was a way to sidestep the extrapolation of the the gridnodes beyond the datapoint boundary. It should be added in matlab next release!

[Darren Paget]

Thanks heaps for the function, it saved me a lot of work. Surface fitting has come a long way! absolutely awesome fit to my data. and really easy to use.

[Henrique Mendes]

It is not working! Matlab says:
"??? Input argument "x" is undefined.
Error in ==> gridfit at 325
x=x(:);"
I tried to change x(:) to x(:,1), but the problem kepps the same. My Matlab is version 7. It had to be more recent?

[John D'Errico]

Henrique - My guess is that you have called gridfit with no arguments. This code is not a gui or command. See the demos for examples of the use of this code. It truly does work, at least when called with data. John

[Gassilloud rémy]

How to use gridfit with polar coordinates ?
[th,r] = meshgrid((0:5:360)*pi/180,0:5:300);
[X,Y] = pol2cart(th,r);
Z=griddata(x,y,z,X,Y);

Thank you,
Rémy

[Tim Marvin]

 

[V P]

What is wrong with calling GRIDFIT? Needs improvement.

Z=rand(100);
a=1:100;
[X,Y]=meshgrid(a,a);
good=(X-50).^2+(Y-50).^2>100;
x=X(good);
y=Y(good);
z=Z(good);
bad=~good;
xi=X(bad);
yi=Y(bad);
zi=gridfit(x,y,z,xi,yi);
??? Error using ==> gridfit at 404
xnodes and ynodes must be monotone increasing

[John D'Errico]

The problem that VP has is he failed to read the help, or apparently bother to understand what gridfit does. I'd recommend reading the help, and perhaps looking at the examples when you don't understand what you are using. I do understand that real men don't need no steenkin help. Perhaps this is the approach that VP has followed.

Gridfit does not interpolate a surface at some random scattered list of points, as would griddata.

Gridfit generates a surface from scattered data on a complete, regular lattice of points. The node arguments are what he might have passed into meshgrid. Again, READ THE HELP.

[V P]

1. In strict accordance with the help, xi and yi WERE x-nodes and y-nodes of interest, generated by MESHGRID. If this kind of input is not accepted by GRIDFIT, this means that GRIDFIT must be improved.
2. This is also necessary for comparison with GRIDDATA, which has no problems with my example.
3. If these arguments are still insufficient, and you insist that I am incorrect, then you have to improve the description of xi and yi in your help and define them as VALID INPUT arguments for MESHGRID. However, this will again mean an improvement.
4. If this is not convincing for you, look at the example below. If it is not a strong indication of necessary improvements, then I have to follow the principle of "haende khokh" - this is much better than the cold war around an empty egg.

Z=rand(100);
a=(1:100)';
[X,Y]=meshgrid(a,a);
z=gridfit(X(:),Y(:),Z(:),a,a);
norm(z(:)-Z(:))
ans =
    4.9225

[Alex Gardner]

 

[Carlos Adrián Vargas Aguilera]

Hi John:
  
I've been using your incredible program for some years and in a wide range of applications. My question is, do you have written any kind of documentation in a scientific journal, for example, for any citation and also to know a little bit more about this super "surface fitting". Any reference can help too.
  
BTW: VP uses the more silly test to "compare" GRIDDATA and GRIDFIT, and, as John says, he doesn't bother to read that by default the former is 'linear' and the latter uses 'smooth' to 1. But, any way...

[Anmol Agrawal]

Excellent utility!! A default setting did everything I needed with minimal error.

[gai cs]

it is very good toolbox ?thanks

[Andrey Matseevskiy]

You wrote "For example, my computer took 1020 seconds to solve a 500x500 problem." What an old comp do you use? I need 5 second to solve 500x500 problem. Download my prog from http://www.smartfills.com/Html/2D.zip . Compare with your own method. My own is in fact interpolation, but could be modified for approximation without much problems.
Andrey, Kamchatka

[John D'Errico]

The review by Andrey is a bit of self serving braggadocio, since he does not actually offer Matlab code on the file exchange that replicates in any form what gridfit does. By the way, the link he gives is for a .exe file - BEWARE.

His time comparisons are also meaningless of course, since the times referred to are for a now wildly obsolete computer, and comparing an interpolation code to this class of modeling code makes no sense anyway.

Finally, converting an interpolation code to a noise reducing surface modeling code is not at all trivial. I welcome Andrey to do so and provide MATLAB code for the purpose, if that is truly so easy.

[Cedric Testaz]

 

[Pauline Wong]

Versatile smoothing function.

[Cristina Velasco]

 

[Alexis]

Good job John, beautiful tool! Very parameterizable and helps a lot in giving you a feel of how your data is behaving.

Does anyone know how I can anchor the surface to any given point? I'd like to force the surface to pass through a given set of (x,y,z) anchor points, something like establishing an "infinite" weight to these keypoints, while not altering the weight of the other "normal" points...

[Jveer]

wow! looks awesome! was wondering if there was any way to fit a surface around a 3D scatter? maybe this could be used on the data of for e.g x,y data corresponding to max(z) and then repeating for max(x) and max(y)?

[Michael Jordan]

 

[J Ú]

 

[Bryan Raines]

 

[Michael Dupin]

 

[Karl]

This is spectacular. I hope the Mathworks is paying you some sort of royalty for your efforts!

[Chris Sherwood]

This routine is almost magic...it is the only routine I have found that finds surfaces I know are there amidst very noisy elevation data. Thanks for contributing it.

[Angel Atanasov]

Thank you mr. D'Errico, for making my life easier.You piece of code is really good.

[Daniel Arevalo]

Hello John, the function is really nice but I'm wondering if it would be possible to use it with non-rectangular domains where the nodes are just specified with (xn,yn) pairs n=1,..., N
Thanks for any help !

[Yaroslav Bulatov]

What is the advantage of using gridfit instead of griddata followed by smoothing?

[John D'Errico]

Use of an interpolant followed by a smooth is a poor second choice, for several reasons.

Gridfit finds the surface that is as smooth as possible, that is consistent with the data. Smoothing a interpolated surface after the fact does not ensure that the result is consistent with the data. When you do a posterior smooth of the surface, the act of smoothing is now disconnected from the data.

Next, if you use griddata to interpolate a surface in advance, you will only get a result that lies within the convex hull of the data. Griddata will not extrapolate unless you use the v4 option, and that option is VERY slow for any significant number of points. Gridfit can extrapolate using several methods, depending upon your goals.

Extrapolation is an important capability of gridfit. But, extrapolation can come in many different forms. For example, consider a data set which is just slightly concave along one edge of the data. The published demo has a good example of this. See the second example, fitting a trigonometric surface. Along the edges of that data, see that the griddata interpolant generates long, thin triangles. Long thin triangles are terrible for interpolation, so what happens is you see strange interpolation artifacts along the edges.

Gridfit allows you to have replicates in your data, treating them properly in a least squares sense to generate the surface. Try out griddata with replicate data points. Even near replicate points can introduce nasty artifacts in the interpolant. Worse, the delaunay triangulation used in griddata will often have problems if you have sets of collinear data points. This is no problem at all for gridfit.

Finally, you can easily control the extent of smoothing done by gridfit.

In short, griddata has its purposes. There are general circumstances when I recommend griddata. But I would never recommend the use of griddata to be then followed by a smoothing operator.

[Michael Adsetts Edberg Hansen]

 

[Ulrich]

 

[Chamane]

excellent

[Christine Acou]

 

[Félix Martin]

 

[Panagiotis]

Is it possible to totally avoid any extrapolation?

[John D'Errico]

In answer to the question of extrapolation, there are several issues here.

Gridfit tries to extrapolate gracefully, that is, it will extrapolate as smoothly and linearly as possible. This is the default mode. It is also possible to specify the 'springs' method for gridfit. This method tries to extrapolate minimally, while still generating solutions over the entire domain of interest. Think of that method as extrapolating as a constant, where it can do so, at least as a continuous, smooth function.

Regardless, these methods still do extrapolation. If you absolutely must avoid any extrapolation, leaving the entire domain outside of the convex hull empty, then an option is to use griddata. Griddata is of course an interpolant. It will do no smoothing as gridfit allows you to control. But griddata will interpolate inside the convex hull, and leave points outside as NaN.

The last possibility is if you need to do smoothing, but only inside a given, arbitrary domain, convex or not. This requires a new tool that is like gridfit, one that I have written but not posted on the FEX yet, but will do so.

[Panagiotis]

John thank you for you prompt reply.

I really need a tool that only interpolates in a given arbitrary domain, as I need to create contour lines in it, which wouldn't make sense out of the domain boundaries.

Could you please post a message here when you post your new tool?
many thanks!

[Jeff Evans]

 

[Matteo Niccoli]

Excellent work! I needed a sound extrapolation algorithm to pad some gravity data to use with both vertical continuation and source depth inversion. It worked very well for me. Might I add this has one of the best and most easily understood helps I've seen for Matlab code, even for someone like me whose first language is not English.
SOmeone asked before about journal reference and I did not see your reply. I'd like to know as well if there's anything published we can reference. Thank you.

[John D'Errico]

Perhaps I should add this...

If you really desire something written, many years ago I was granted a pair of US patents on a similar idea, there applied to higher dimensional problems in color modeling. Gridfit only works in 2-dimensions of course.

US Patent # 4992861, 4941039

At least for me, actually reading a patent is about as miserable a pastime as I can imagine.

[John D'Errico]

Alternatively, a simple citation style for a FEX submission is found here, if all you wish is a citation:

http://matlabwiki.mathworks.com/Citing_Files_from_the_File_Exchange

[Janine Bijsterbosch]

Hi there! Is there a gridfit3d available at all? I have a 3d array of electric field data (369x551x325) on a irregular but highly collinear grid and griddata3 therefore does not work. As I understand, the gridfit code would be perfect for this if it could deal with 3d data. Thanks!

[John D'Errico]

I've responded via direct e-mail, but the gist of it is that I have plans to introduce a gridfitn one day, but it will take some time to write, and I have many things to write on that same list.

[Jeff]

 

[Adam Chapman]

This is brilliant, but a capability to vary the degree of smoothing in x and y exclusively would be very useful

[John D'Errico]

I like Adam's idea, and this is easily enough put into the code.

[Adam Chapman]

the new setting "'smoothness',[xsmooth ysmooth]" is a Godsend. Many Thanks

[xfactor y]

 

[Ben]

Thanks John for your code. I enjoy it very much!

[Alexis]

+1 vote for gridfitn ;)

[Lars Robert]

excellent job on meteorological balloon data!

[Toni]

Hello!
I have an 512x512 array z where ~11000 points are non zero and the rest must be nans! Why is the following code giving me an array zgd filled with only nans? How does this program deal with nans?

xi = linspace(0,1,512);
[xg,yg]=meshgrid(xi,xi);
x = 1:512;
y = 1:512;
zgd = griddata(x,y,z,xg,yg);

[John D'Errico]

Argh.

Gridfit is set up to ignore any nan data. So your question is meaningless and irrelevant, because gridfit in fact ignores the nan data. See the following code fragment, taken directly from the code...

% also drop any NaN data
x=x(:);
y=y(:);
z=z(:);
k = isnan(x) | isnan(y) | isnan(z);
if any(k)
  x(k)=[];
  y(k)=[];
  z(k)=[];
end

So giving this code a rating of only 4 stars because a different code (griddata) did not work for you seems a bit silly.

If your goal is to fill in the nan elements, while leaving the actual data alone, then use inpaint_nans to interpolate. Use the right tool to solve your problem.

If your goal is to do smoothing, while also interpolating the empty (nan) data, then use gridfit. It will solve your problem, but only if you call gridfit and not griddata.

Finally, if your problem is that griddata is not working as you wish or that you don't understand griddata, then why in the name of god and little green apples, why are you asking it here, in a place intended for comments on gridfit?

[Toni]

Sorry!
Still stumbling about these issues.
Thanks for the right directions!

[Michael Brown]

Just wanted to extend my gratitude that you made this available. The documentation is excellent, it works as advertised and saved me several days of setting up and debugging a regularization scheme for fitting an appropriately smoothed surface through some unequally spaced thermodynamic data that have error and do not extend to all corners of the needed P-T grid. It was a relief to be able to focus on the science rather than on debugging the scheme to process the data.

[Lennert]

Hi there! I have a question in relation to this tool. Suppose you have already fit a surface through the data. Is there a command to extract the formula/function of that surface? Thanks!

[John D'Errico]

There is no functional form for the surface, at least unless you are willing to accept a tool like interp2 to interpolate it at any point. The best that you can do is to view that surface as a low order spline in two dimensions, but as a spline, all you have are a large number of tiny linear segments all neatly connected together. (This applies to the default method, which uses a triangulation of the regular lattice. The alternative is the bilinear interpolant, which in effect is not even truly linear, but an oddly and mildly piecewise quadratic form.)

The virtue of this approach is you can fit any set of data that follows the general form of a function z = f(x,y). You need not formulate a model as other modeling tools force you to do. But you can't get a model out of this either.

If you really want to do empirical modeling, you need to make the effort of posing a model, of choosing some form that will represent your surface. That model may be polynomial (as my polyfitn would allow you to fit) or it may be nonlinear. Then you need to use a tool that can fit that model to your data. There are many such tools on the file exchange to serve that purpose, or you can use the curve fitting toolbox.

[Jason]

I have scattered 3d data of a membrane subject to tensile loads at the corners resulting in transverse displacement. I'm trying to use gridfit to look at a the 3d surface representing the data. The problem is around the edges the representation gridfit and surf produce is very inaccurate. In fact, at the corners the resulting plot shows negative values when there are no negative z displacements in the data! Any ideas? Thanks.

[John D'Errico]

Jason - This is a problem of extrapolation, something often difficult to do intelligently. My standard response is to quote Mark Twain:

“In the space of one hundred and seventy six years the Lower Mississippi has shortened itself two hundred and forty-two miles. That is an average of a trifle over a mile and a third per year. Therefore, any calm person, who is not blind or idiotic, can see that in the Old Oölitic Silurian Period, just a million years ago next November, the Lower Mississippi was upwards of one million three hundred thousand miles long, and stuck out over the Gulf of Mexico like a fishing-pole. And by the same token any person can see that seven hundred and forty-two years from now the Lower Mississippi will be only a mile and three-quarters long, and Cairo [Illinois] and New Orleans will have joined their streets together and be plodding comfortably along under a single mayor and a mutual board of aldermen. There is something fascinating about science. One gets such wholesale returns of conjecture out of such a trifling investment of fact.”

"Life on the Mississippi", Mark Twain, 1884

The point is, extrapolation is difficult, especially if all you use is a tool that tries to take data that it sees and predict into the unknown. How should matlab know that zero is a special number? How should gridfit know that it is ok to predict smoothly beyond the extent of your data, yet stop at zero?

Having said all of that, there are several options open to you, both of which might work nicely. First, you can use the 'springs' method in gridfit. It tries to prevent extrapolation beyond your data. This is sometimes the proper solution. Try it, and you might be happy, or not. This I cannot say.

A second option is a transformation. Very often when a system has a property that it cannot go less than zero, you are working in the wrong "space". The trick here is to use a transformation to fix it. Here, I might log your data. Now, use gridfit to model the surface, allowing it to extrapolate smoothly as it wishes with the default options. Now, exponentiate the result. In effect, you have done an interpolation in the log domain, but then transformed back. Do you care if the logs went negative? Of course not. exp is a function that is positive for all real inputs. This trick often works very nicely. In effect, it is just recognizing that the interpolation is best done in the proper domain, one where your system is truly more additive.

[Mikkel Lund]

Very nice.. Just what I have been looking for

[jose tapia]

Hi,
I'm having the following problem. I have a 2d table that contains data from some measurements. How can I do the same kind of Extrapolation that is possible in SIMULINK 2-D table lookup using interpolation-extrapolation lookup method, but in Matlab. As I figured out 'griddata' and 'interp2' can not do the job for me.
This is the dimension of my data:
x <1x37 double>
y <1x28 double>
z <28x37 double>

Thanks in advance.

[John D'Errico]

You have asked this question in multiple places. I have no idea what simulink does, so I cannot answer your question, nor can I even be sure what it is you wish to do.

Matt already answered (as an answer) that if you just desire to interpolate the function at a specific location, then interpne (from the file exchange) will solve that problem. I believe that is the correct answer from what little I know about your problem.

If you wish to create a new surface that extends further out, gridfit can do that. Just supply the new coordinates that include the old x and y, but will create a new grid that extends out as far as you wish.

Remember that any extrapolation is a risky business, trying to predict the unknown. It is RARELY accurate or even intelligent, and is only as good as the data it is built on. If you are trying to extrapolate a long distance out, don't be upset at the result.

[jose tapia]

Thank's for your answer John D'Errico.

[Ignacio]

Really good

One way of using gridfit just for interpolation is :
- Use gridfit in the whole grid zifit=gridfit(X,Y,Z,xi,yi)
- Use griddata in the same grid
zidata=griddata(X,Y,Z,xi,yi)
- Use isnan to know which elements are NaN in zidata, and apply it to zifit

At least for me it worked
- Use a isnan(griddata) to know which points are inside the original data, and apply it to the points created with gridfit

[Januka]

Great piece of work John. No hassle and easy to use. Thanks!

[John Reinert Nash]

 

[John Reinert Nash]

Excellent tool that much improved my visualization of some photographic exposure data. Thanks, John, great to use one of your tools (again).

[Melissa]

This is an amazing piece of work John. I bow to your superior programming skills and knowledge. I have a question though, is it possible for me to evaluate a point or a series of points on the surface?

[John D'Errico]

Melissa - gridfit merely generates a function defined at a set of discrete node points, on a regular lattice. It does not deliver an interpolated prediction at any point. But you can easily enough get that prediction from the surface produced from gridfit.

There are three interpolation methods that are allowed in gridfit. Those methods define how gridfit treated your data when it built the surface itself. Nothing requires that you use the same method to interpolate that gridfit used when it built the surface though.

The nice thing is, we already have interp2 to do 2-dimensional interpolation once the surface is defined on a regular lattice, with several good methods already there. The case of 'linear' for interp2 is equivalent to that employed in gridfit for the 'bilinear' case. It turns out that this is also a method used by tools like photoshop for image interpolation.

You can also use the other methods of interp2 though - splines for example. Or, if you prefer a simplicial interpolant, where each square of the lattice is split into a pair of triangles, then a truly linear interpolation is done across those triangles, you can find my interpns on the file exchange.

http://www.mathworks.com/matlabcentral/fileexchange/30932

So yes, you can easily do interpolation at any arbitrary point that lies inside the lattice as defined by gridfit.

[stefano pasquali]

hi John.
i'm pretty new in matlab. so maybe i'm going to do a stupid question.

i have my data already in a matrix form.
tipically 200x200.

i cannot use your function because it expect 3 vectors and nodes ...

can you help me ?

basically i have this surface, results of observation (along x and y there are categorical variables) and i'm expecting a monotonic decreasing surface ... i got it but with some outliers ... i think your function can help me.

any suggestion ?

[Qiang]

hi John,
I got error information now, but not last year. Seemingly it is because the code is not consistent with the new Matlab version (2011a). Maybe I am wrong. Do you have ideas about this? The error is as following:

??? Error using ==> mldivide
MLDIVIDE is not supported for one sparse input and one single input.

Error in ==> gridfit at 616
        zgrid = reshape(A\rhs,ny,nx);

Thanks

[Yang]

Hi John, I have a set of 2d scattering data, and I want to have interpolation values at a regular (1:N,1:N) lattice. What parameters should I use to minimize the smoothing? (I need to iterate the interpolation for 50+ times, setting "smoothness"=1 for 50 times makes the resulting image very blur).
Thanks!

[John D'Errico]

Smaller values of the smoothing parameter give less smoothing, although zero as a value may result in a singularity.

[Yang]

Hi John,
I tried some smaller parameters (say 0.1 for smoothness), but then I get wave-like artifact pattern in my image.
What I am looking for is a combination: griddata in the inner part (so no smearing at all), and gridfit at the edge (where griddata produces NaN). Hence smearing is minimized. Can I achieve that by tuning parameters? Thank you very much!

[M. C Ertem]

Thank you!!! Saved the day. (Actually, night.)

[Christian]

I'm trying to use gridfit to produce a contour map (functions: contour and contourf). MATLAB seems to interpret the griddata output okay for the surf() function, but fails with contour(). Anyone else run into this problem? I think it has something to do with grid orientation.

[John D'Errico]

How can anybody guess what you are doing wrong? How exactly does it fail? What are you doing exactly? What is your data? And why are you saying that griddata output is a problem, in a question about gridfit?

[Christian]

I'm sorry - that was a typo: I meant gridfit when I typed griddata. I since figured it out it was a MATLAB environment issue, unrelated to gridfit. Contour plots just weren't showing up in the display. Restart solved the problem. I would edit or remove my original post, but this forum doesn't allow for it. Thanks for the quick reply.

[Darin]

Thanks for the great tool. Trying to find a way around "pull ins" when the regularizer has a significant low spatial frequency bias over extended areas... for example, around the top of a higly resolved gaussian peak with 'gradient'. Using a higher order regularizer doesn't help (I've implemented a 4th order gradient): that just gives a surface that's highly "puckered" around indiviual sample errors. Any ideas?

[Mark]

 

[Mark]

Very nice, currently using this to perform a calibration between an imager and a gimbaled laser pointer. It's not quite optimized to get the best expected performance over the full FOV, but it's close. Toward the edges, it works much better than other routines I've tried. Is there a way to see how good the fit is?

[John D'Errico]

This will definitely have much better performance around the edges compared to a delaunay based tool. They often show serious interpolation artifacts around the edges.

In terms of how good is the fit, I'm surprised that I did not return an array of residuals or predicted values, but this is easy enough to build. Use interp2 to compute predicted values, since what comes out of gridfit is now a nice regular one that interp2 will handle. Once you have the predictions, residuals are easy enough to get by subtraction, or you could compute a standard error, etc.

[Peter]

GridFit is pretty good, but I have a few questions.

1) The interpolation seems to be triangular, which causes some symmetry problems. Below is an example program to illustrate this. Does anybody know if there is an advantage to using triangular interpolation (three points) instead of rectangular interpolation (four points)?

2) If I change the grid spacing, I also have to adjust the smoothness parameter or GridFit gives different results. This appears to be related to the way GridFit sets up its derivative equations. It mixes horizontal and vertical second derivatives in the same equation. For example, if the residual in the horizontal direction is +0.1 and the residual in the vertical direction is -0.1, the regularizer wrongly deems that point a "perfect fit" because the derivatives are added together. Would it be better to separate the derivatives into separate equations (separate rows in Areg)? If so, is it possible to produce the same results (same smoothness or curve profile) for a single "smoothness" parameter, no matter what grid spacing is used? I think that would be really helpful and convenient to GridFit users so that we wouldn't have to fiddle with the smoothness number manually.

% Choose four points at each of four corners and one point in the center.
% Make sure everything is symmetric around the center.
SourceData = [
  0, 0, 0
  0, 1, 0
  1, 0, 0
  1, 1, 0
  0.5, 0.5, 1
  ];

% Create a 4x4 grid.
xNodes = linspace(0, 1, 4);
yNodes = xNodes;

% Run GridFit. The smoothness parameter has little effect on the result
% for this demonstration.
[zGrid, xGrid, yGrid] = gridfit(SourceData(:, 1), SourceData(:, 2), SourceData(:, 3), xNodes, yNodes, 'smoothness', 0.1);

% Display the surface so that we can see the asymmetry.
surf(xGrid, yGrid, zGrid);
xlabel('x');
ylabel('y');
zlabel('z');

% Calculate the asymmetry.
disp(['GridFit asymmetry is ', num2str(zGrid(2, 2) - zGrid(3, 2)), '.']);

[John D'Errico]

By default, it uses triangular interpolation. If you read the help, you will see that a bilinear (quad based) interpolation method is an option. Is one better than another? In general, I'd claim that if you can see the difference, then your grid is far too coarse! In any event, HAD you actually read the help before posting, in your example you might have tried adding one more option to the call.

[zGrid, xGrid, yGrid] = gridfit(SourceData(:, 1), SourceData(:, 2), SourceData(:, 3), xNodes, yNodes, 'smoothness', 0.1,'interp','bi');

Your second comment is about the smoothing parameter. Note that gridfit allows a common scheme for the smoothness penalty function - that the Laplacian is biased towards zero. This is the sum of the second partials, and it is arguably the logical choice for SOME physical surfaces. HOWEVER, note that the default is NOT that method. In fact, the default is a method that DOES uncouple those second partials!!!!!!!!!!!!!!!! Again, read the help, rather than assuming you know the answer and then asking a question based on that wrong assumption.

As for the optimal smoothness parameter changing based on the grid spacing, I have found that for virtually any method of choosing a smoothness parameter, I can also come up with a case where it will not be the best. Gridfit uses a default smoothing parameter that is "reasonable" for many problems out of the box. Is it optimal? Probably not. Any degree of true optimality might involve algorithms that would be slower to run, and algorithmic speed is desirable here. Should there be an adaptive option in gridfit, that would allow the user to set it and walk away, knowing that it will be slower to run, but it will give an always optimal result, for all users? Truly good adaptive methods that will never fail are also terribly difficult to write. Look at the numerical integration tools in MATLAB (quad, quadl, quadgk, etc.) I can easily make those tools fail by a careful choice of integrand and interval.

All that I can suggest is to use the facilities built into gridfit as a computational tool. If it produces by default not exactly what you want, then I have provided a few knobs you can adjust. Only you know if you like the surface you create with a tool like gridfit. But at least the knobs are there to turn to generally produce something that will make you happy.

[Yan]

Hi, this is an awesome program! But it seems to be very memory consuming. I can solve for Zgrid points with a rather small dimension only (e.g. 300x300). Griddata seems to be able solver for a larger matrix faster, but the interpolation is not nearly as good. Is there anyways to obtain higher resolution?

[John D'Errico]

What can I say? You can have a good result, or a faster one. Pick only one.

The fact is, gridfit is actually blazingly fast given what it does, which in this case requires the solution of a linear system with 90,000 unknowns. A quick test shows that took only about 2 seconds on my system. Perhaps you want me to do magic? Even a wave of a magic wand takes about 2 seconds.

[Tomas]

Great work!