Model 2d surfaces from scattered data
Editor's Note: Popular File 2008
This file was a File Exchange Select file.
Select files are submissions that have been peerreviewed and approved as meeting a high standard of utility and quality.
Those wishing to model a surface from data in the
form of z(x,y) from scattered or semiscattered
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.
1.1  Added a live script 

1.1  Added a live script 

1.1  Allow unequal smoothing parameters in x and y directions 

Release 2.0: Adds the option of tiling to solve much larger problems 

Allows the user to specify a number of nodes from

Inspired: RegularizeData3D, smoothsurf( x,y,z,xn,yn,varargin ), jasonnicholson/regularizeNd, Bi dimensional Emperical Mode Decomposition (BEMD)
Mahdi S. Hosseini (view profile)
A relevant toolbox for numerical differentiation called MaxPol has been recently published and made available here
https://www.mathworks.com/matlabcentral/fileexchange/63294maxpolsmoothinganddifferentiationpackage
MaxPol provides a framework to design variety of numerical differentiation kernels with properties like:
(1) Cutoff (lowpass) design with no sidelob artifacts (for noiserobust case)
(2) Arbitrary order of differentiation
(3) Arbitrary polynomial accuracy
(4) Derivative matrix design
(5) 2D Derivative Kernels with Steering moments
(6) Intuitive examples in Signal and Image processing
xiaojie chen (view profile)
Howerver , I don't understand the difference between this one and "scatteredinterpolant" in the function hub of maltab
Gert Kruger (view profile)
Jason Nicholson (view profile)
If you like this function, you should try regularizeNd. It is nD version of gridfit with some bug fixes and improvements.
https://www.mathworks.com/matlabcentral/fileexchange/61436jasonnicholsonregularizend
There is a bug related to the number of 2nd derivative equations added to the fidelity equations, A. The effect is that have added equations to the system that look like 0*x = 0 when trying to solve for x. Because it is a least squares problem, this equations doesn't destroy the solution and ends up only causing minor issues but it is a bug.
Also, gridfit doesn't scale the 2nd derivative equations. This is important when the units of 1 dimension are different than the other or the scales are very different. The smoothness parameter needed in one dimension should be similar in value but they are not in gridfit when dimensional scaling is different.
Geoffrey Vincent (view profile)
Ernst Niederleithinger (view profile)
Arun Krishnadas (view profile)
John! This was very helpful. My huge vectors imported from abaqus results were struggling to get fit through griddata. gridfit helped solve the issue.
Sebastian Bannwarth (view profile)
Dear John! What a great tool. It helps me a lot fit a grided surface to my scattered data. I wonder if there is any opportunity to create a sfitobject with gridfit? What i actually want to do is differenting the fitted surface at some spot and later integrate again. So is it somehow possible to get the surfave as an sfit object?
Guojin Feng (view profile)
Sven (view profile)
John, as usual a wonderfully useful utility applicable to so many problems.
I thought I might ask  can you see any way that periodic boundary conditions may be implemented into the workings of gridfit?
I have encountered this problem on a few occasions and so far I've just gone with some hacks to approximate periodicity. For example, if in my data I expect periodic output along the x direction (such that the first and last columns of gridfit's output surface should ideally match each other in position and slope), I've just taken the regular gridfit output, padded that output with mirrored copies of the opposite few columns, then applied some gaussianlike smoothing to those edge columns before removing my padded additions.
As I said, it's always been a hack and is never perfect  but something like this has often been "good enough" for my purposes at the time.
So, can you comment on the feasibility of adding constraints within gridfit to match the position/slope of the first and last columns (or rows)?
B Verhaar (view profile)
Dear John,
Thanks for this wonderful function. I make use of it a lot.
It appears there is a 'hidden' input argument, called mask. What is the purpose of it?
Kind regards, Boudewijn Verhaar
Jianguang Chen (view profile)
Excellent work!
Binu (view profile)
Somayeh Hesabi (view profile)
pengfa zhou (view profile)
hi,i'm a learner of matlab, anyone can tell me where i should paste this file for use?thanks million time!
j b (view profile)
talyn dong (view profile)
Natalylun (view profile)
Serge (view profile)
RegularizeData3D is a fork by Jamal where he has:
1) fixed normalisation issues, such that chaning grid resolution does not change the shape of the fit.
2) added bicubic interpolation which can make a small difference in some cases.
Sadly, neither version works on 2D data. Some hacking is required to create a fake 3rd dimension with something like [0 1 0 1...]
http://au.mathworks.com/matlabcentral/fileexchange/46223regularizedata3d
rajyalakshmi k (view profile)
hi,
i have xyz values i,e .mat file, i need to reconstruct surface with them i am using this code but getting errors can anyone suggest me the procedure how to use this code properly.
Thank you
Ehsan golk (view profile)
it is not working for nonmonotone increasing vectors
Behailu Shikur (view profile)
Stephanie (view profile)
Hi!
I'm just come across this and it looks really great!
I was wondering if this also works for segmented surfaces, like two surfaces intersecting each other at one point.
Would that be possible?
I could not figure it out how to use it in that case.
Thanks for any help!
Stephanie
Lyra Maxwell (view profile)
Julia Howe (view profile)
This is a great tool  much more efficient than griddata.
I was wondering if there is a way to extract the formula of the fitted surface using this code?
Thanks,
Julia
Zhang Lidong (view profile)
Very Good
Rose Ann Haft (view profile)
Does this create a formula of the fit?
Sathish Sanjeevi (view profile)
yousaf obaid (view profile)
arnold (view profile)
Hi John,
after using this very robust function for years I finally stumbled upon a question I can't answer myself and that is how to use weighted points. Since this is not implemented I wonder if anyone here could point me to an alternative.
I see Matthias was working on something in march?
kind regards
Arnold
Eugenio (view profile)
I use it often, thanks for your job
INDRANIL SAAKI (view profile)
Pietro (view profile)
Rob Campbell (view profile)
daniel perez (view profile)
Hi,
When using the gradient regulizer, why does
the doc file says gridfit attempts to force the ''first'' partial deratives to be equal around a node? When I look at the code, it seems to me that is forcing the second partial deratives to be equal to zero around every cell. Not exactaly like the laplacien method indeed, but still dealing with the second derivative.
code where I see a second derivative under the gradient regulazer:
2./(dy1.*(dy1+dy2)),2./(dy1.*dy2), 2./(dy2.*(dy1+dy2)
Jan Mervart (view profile)
daniel perez (view profile)
Amazing work!
Henry1993x (view profile)
Wrichik Basu (view profile)
Matthias (view profile)
I tried implementing a weighted version myself. I did it by replacing the Normal equation in the "normal" solver (gridfit line 631) by the weighted equivalent
zgrid = reshape((A'*W*A)\(A'*W*rhs),ny,nx);
where W is a diagonal matrix with the point weights on the first nPoints entries of the diagonals and ones on the rest. The weights are normalized so that they sum to nPoints. Is that a reasonable approach?
Matthias (view profile)
Thank you very much for making this great function available!
I'd like to be able to apply different weights to different data points and am wondering if there is a rigorous way to do this. I guess I could create copies of data points according to their weights but that seems like an inelegant solution.
Thanks!
samar (view profile)
Can someone help with this error?
Undefined function 'sparse' for input arguments of type 'single'. Thanks.
Binu (view profile)
Matthew (view profile)
Great submission as always John. Thanks.
Already asked below, but can you add bicubic interpolation to this submission?
Please :)
Giorgio (view profile)
Erika (view profile)
PONYO (view profile)
bigtoelee (view profile)
hey thanks very much for the code.
i have a question:
my 2d data are formatted in matrix already, do i have to convert them into column to use the function?
MaLe92 (view profile)
Okay, thank you.
John D'Errico (view profile)
Think of the result from gridfit as essentially a spline. This is really what it is, a piecewise function. Piecewise functions have no simple representation of the form f(x,y). At best, you could write that evaluation in the form of a call to interp2, because what gridfit returns is your function over a rectangular array of points.
Finding some general function of the form f(x,y), in a nice form that you could write down for a paper, etc., is a difficult problem to do (automatically) by computer. Such a function would in general be a nonlinear function of x and y, and it would involve parameters that must be estimated. To do this, you first need to decide on a nonlinear model, then use a nonlinear estimation tool to determine the parameters. None of that is in gridfit, nor can it be.
MaLe92 (view profile)
Hey, :) Thank you very much for this code.
I have another question:
Is it somehow possible to determine the generated surface as a function z=f(x,y)?
I can't find information (or i don't understand the information ;) ).
John D'Errico (view profile)
No other implementations that I know of, at least that are available for free. I have heard of others asking if it was ok to write it in some other language, but I've not seen any implementations pop up. This does not say nobody has written one, only that they may have chosen to keep it out of the public domain. Given that it may take some effort to get it working, I would not blame someone who wrote such a code but kept it to themselves. For example, I wrote an ndimensional version of this long ago, but that code belongs to my past employer.
Jarrod Collins (view profile)
This code is wonderful.I was wondering, is there a C++ or VTK implementation to this?
John D'Errico (view profile)
Think of the surface produced by gridfit as a moderately flexible plate, where we can control the overall flexibility. Now, attach springs that attach to the plate, connecting the plate to each data point. Springs have the property that the potential energy stored in the spring is proportional to the square of the extension, so a least squares kind of thing. Double the extension, and you multiply by 4 the energy stored in the spring. The springs are designed to store energy only in the zdirection.
The basic idea to the smoothing parameter is that we can adjust the flexibility of the plate relative to the spring constant for those connecting springs. So if we make the springs stronger, then they will distort the shape of the plate more. Make the springs very weak, then the plate will approach a planar surface.
GdogCdog (view profile)
Excellent code  speedy and works straight from the box! What more could you ask for.
One question, Could you explain a little more about how the smoothing parameter works. I have tweaked to get the desired result, but want to understand, and be able to explain to others a little more about how this works.
Thanks!
Junghwan (view profile)
Thanks for the great code.
I have the same request as Jia Le Ngai.
It expands over the data and fills the whole rectangular area. Is there a way to stop it extending?
John Porter (view profile)
I have used this program for years and it works well! Has anyone tried converting gridfit to C using Matlab Coder ?
John D'Errico (view profile)
Gridfit is not an interpolation tool. If you want that, use tools like griddata, triscatteredInterp, scatteredInterpolant, etc.
Jia Le Ngai (view profile)
Thank you John for the code. However Is there a way to make the just to interpolate the coordinate that I input without the 'extend' property?
Thank you,
Jia Le
Jia Le Ngai (view profile)
John Chen (view profile)
The code is great tool and the discussion is very helpful. I have run demos, the results are very impressive.
I notice you have another toolbox which able to fit surface on spherical coordinates. Could you please send me a copy of it?
Thank you very much!!
Neil (view profile)
Thanks so much. I am very appreciative.
Neil
John D'Errico (view profile)
Whatever value of z0 that seems to work should be fine there.
You can think of it as a truncated Taylor expansion if you like, which it essentially is. I simply chose a straight line that matched the slope and value of the log function at the break point, to create a continuous, differentiable function that will extrapolate linearly. So sort of a spline too.
Neil (view profile)
Thank you, John; your suggestion worked. I chose Z0 to be 0.5% of max(Z) of 2E12 in the data I'm fitting. Now, using gridfit to fit the conditionally log transformed data yields a positive surface which does not overshoot.
A question: how did you arrive at the second part of the transform, i.e. ((z  z0)/z0 + log(z0)).*(z>z0)? It works but it's not immediately obvious to me why. (Some sort of Taylor expansion of the log?)
Thanks,
Neil
John D'Errico (view profile)
Neil,
Of course, my response should be  picky, picky, picky!
This is a consequence of a log transform I guess. You could use the 'springs' method as a regularizer to avoid that, but it might be too extreme.
So an alternative is to use a piecewise transformation. Thus, choose some breakpoint, above which the transformation will be linear, and below which, it will behave as a log transform. I'll call that value z0. So we can write the transform function W(z) as
W = @(z,z0) log(z).*(z<=z0) + ((z  z0)/z0 + log(z0)).*(z>z0);
You can plot it, picking perhaps z0=3 as a breakpoint for an example.
ezplot(@(z) W(z,3))
Of course, the transformation is invertible. This should work:
Winv = @(w,z0) exp(w).*(w<=log(z0)) + z0*(w  log(z0) + 1).*(w>log(z0));
Again, we can plot it, here for z0 = 3 again, as:
ezplot(@(w) Winv(w,3))
Something like that should give you the best of all worlds, and all you need to do is choose some value for z0 that makes sense to you. I might start with the mean of z as a good choice for z0.
Neil (view profile)
Hi John. Thanks for the tip regarding using a log transform. It seems to be a sensible approach and I tried it. The fitted surface is now all positive, as expected. However, it also overshoots: the highest peak (Z) in the fitted surface is now 3 times the largest data point. Previously the fit closely followed the data.
I should add that the data did have 0's in it and I didn't want to completely remove them because they represent boundary conditions for the somewhat sparse data I am fitting to; for example, the spectrum vanishes at a certain energy (e.g. E=0) or at some angles. Instead, as you suggested, I replaced the 0's with a small number; I tried 1E10, 1E5, 1E1, and 1. The overshooting appeared in all cases and is probably unrelated to the presence of 0's or small numbers in the data.
I tried to attach pictures of the fitted surface before and after the log transform, but there doesn't seem to be a way to do that on this forum.
Thanks,
Neil
John D'Errico (view profile)
Neil  the answer is definitely yes, and no. Or maybe it is the other way around.
No, it would be difficult to make gridfit work as you would like to see via internal changes, since that would be a bound constraint on a fairly massively large scale linear algebra problem. So while I could in theory use a tool like lsqnonneg or lsqlin instead of backslash for the solve, it would take forever to terminate on any problem of reasonable size.
Regardless, there is a solution that will work for you, and is quite simple. Typically when there is a bound constraint at zero because of physical issues like this it is because the problem really should be transformed to make it more linear. The logical transformation is a log transform. Thus instead of fitting Z as a function of X and Y, fit the surface in the form of log(Z) as a function of X and Y. Clearly there will no longer be any negativity issues, because to recover the surface you really want, you will simply exponentiate.
If you actually had any data points that were an exact zero, it would be best to drop them completely. Alternatively you could replace them with some small number, but that would cause the problem to be less smooth and might introduce bugs in the surface. Of course, the nonsmooth region will be in areas where the log will be very negative, so after exponentiation it will be essentially zero anyway. So either way will work on those exact zeros.
I hope the log transform idea solves your problem. It is an approach I have used often with good success.
Neil (view profile)
Thank you for making this code available. Is there a way to tweak it to add the additional constraint that the fitted surface should also always be positive? The X,Y,Z data I'm fitting represent spectrum (Z) vs. energy (Y) and angle (X). Therefore, negative Z is unphysical.
Thanks,
Neil
Darin (view profile)
Not only a great tool, the source code and the discussions here are a great education. Thanks, John, for extending the docs and examples
It usually does everything I need. When really offthewall needs have come up, it's been a great base to start from. Clear and well structured code.
John D'Errico (view profile)
Of course, the simple answer for Chad is to put a basic wrapper around gridfit, one that tests the data, and compares it to the nodes in advance. If the nodes were not chosen to contain the data, then fix them so they are. No warning need ever be generated then. Then just have the wrapper code pass all arguments into gridfit.
John D'Errico (view profile)
I (or you) could change the code as Chad has suggested. Personally, I don't think this is a good idea, since gridfit appends new nodes, expanding the grid when it finds this condition. That means the size of the grid as an array is no longer as you might have expected it, so a warning seems merited. And even if the nodes were just moved without telling you, I think this is a bad idea to say nothing.
As it is now, if you really positively don't want a warning, then you can turn that warning off. That is a documented option for the 'extend' parameter. Essentially, as I have it written now, if the code does something to the grid that is not as you expect it, it generates a warning unless you tell it not to generate a warning. I'm just not comfortable with the idea of a wishywashy warning, that may or may not generate a warning when it encounters something surprising.
Chad Greene (view profile)
My xnodes and ynodes are typically ~1 km apart, but if gridfit needs to extend the domain by a nanometer in any direction it'll give a warning message. I recommend adding a clause which will limit warnings to only cases where the domain is extended by more than half of a dx or dy. In the switch params.extend statement there are four possible warning cases. The warning clause I'm suggesting would look like this:
case 'warning'
if xmin<xnodes(1)abs(mean(dx)/2)
warning('GRIDFIT:extend',['xnodes(1) was decreased by: ',num2str(xnodes(1)xmin),', new node = ',num2str(xmin)])
end
xnodes(1) = xmin;
And similarly for the other three warnings the if statement would be
if xmin>xnodes(1)+abs(mean(dx)/2)
if ymin<ynodes(1)abs(mean(dy)/2)
if ymax>ynodes(end)+abs(mean(dy)/2)
Judith (view profile)
This is a fantastic function and it runs perfectly with my 2Ddata. Do you also provide a gridfit3d version? That whould save a lot of time in my postprocessing work. Thanks in advance.
John D'Errico (view profile)
http://blogs.mathworks.com/community/2010/12/13/citingfileexchangesubmissions/
Chaipichit (view profile)
How to cited your code in research paper?
John D'Errico (view profile)
Evgheny  If you seriously need to fit such a surface in spherical coordinates, I might be able to help you out with a separate, unposted toolbox that is capable of fitting a surface in spherical coordinates. You would need to contact me directly.
Evgheny (view profile)
Thanks, function does do its job great!
Does anyone know alternative of this function for closed surface (spherelike).
If I try this function with spherical coordinates, it works fine, but there are problems with the seam (line of joint)
John D'Errico (view profile)
Sorry. I've never heard of a Java implementation. Not that it counts, I've written it 7 times in MATLAB, once each in Fortran and APL, those last were very old though. The difference each time was what I learned from the previous incarnation. Were someone to do a Java implementation perhaps the most important thing is to use sparse linear algebra capabilities for the solve, else it would take forever.
behruz (view profile)
is there any implementation of this in java? that get some points location and value and give a grdi?
Alex (view profile)
Just to add my $0.02 to this whole sphere discussion. By math definition of what is function it means for each x we have only one F(x). If there are more than one F(x) relation between x and F(x) is not a function.
Klont (view profile)
I wish I could remove my previous question, 2 seconds on google pointed me to Matlab's interp2 which I'm using now on the output of your excellent gridfit. Thanks for sharing!
Klont (view profile)
I found that it is very straightforward to fit a nice looking surface to some 55 datapoints using gridfit. I was wondering though, is there a way to get values of the surface that are not on the nodes? I wish to interpolate between my 55 samples using the surface.
Klont (view profile)
John D'Errico (view profile)
For those of you who don't appreciate why Felon's comment is silly, think of it like this. Gridfit fits a surface of the form f(x,y), over a rectangular grid. It does so quite well, as many people have found over the years.
While you may think of the surface of a sphere as a surface, it is not of the form that gridfit can fit. It is multivalued, so for any single (x,y) pair, there will be zero, one, or two values of z that would apply. As well, that "surface" (better to call it a manifold) has derivative singularities, if we were to look at it as a function of x and y. So even a hemisphere will be problematic for this tool.
You would not use gridfit to fit something that is not representable as a function of two variables over a rectangular grid, any more than you would expect it to do numerical integration, numerical optimization, or compute an FFT. Nor would you expect it to cook dinner for you, do your laundry, etc. Use the right tool to solve your problem, but if you try to force the wrong tool to solve a random problem, expect poor results and don't complain about what you get.
John D'Errico (view profile)
Felon 
So you are using a tool that builds a single valued function to fit something that is obviously not. What did you expect? Magic?
Software does what it is programmed to do. It does not magically rewrite itself when you give it a problem of a completely different sort. In fact, I fail to understand why you would downrate a tool for not solving a class of problem it is explicitly not designed to solve.
If you have a closed manifold, like a ball or some other multivalued form, then don't use this tool. I have NEVER claimed it would solve that problem. Instead, you might look into tools like convex hulls, alpha shapes, CRUST, etc. Or, you might choose to convert the problem into spherical coordinates, at which point gridfit would be able to build a viable surface.
Or maybe you just wanted to complain with no good reason.
Highman Field (view profile)
Still it cannot surface a cluster of points that looks like a ball. The result the interpolation is always a surface of singlevalue function.
Highman Field (view profile)
Angelo Tafuni (view profile)
John D'Errico (view profile)
No, gridfit does not explicitly allow you to apply derivative constraints. That does not say it is impossible, only that I did not offer it as an option.
The main reason why not, is it would require a set of linear inequality constraints on the unknowns. For a not uncommon grid of size 100x100, there are 100*100=10000 unknowns to solve for. This is not a problem, since the linear system is a sparse one. However, to solve a sparse linear inequality constrained system, one would need to use LSQLIN, or a solver like it. And the last time I checked, LSQLIN was not set up yet to handle sparse large scale inequality constrained problems. (That may have changed with the most recent release, but I have not checked.) If I made all of the matrices full ones, the solve time would probably be incredibly slow and memory intensive.
So I'm sorry, but gridfit will not handle the problem as is.
If you were willing to build a fairly coarse grid, AND add the constraint system, it would probably be doable in a reasonable time. I don't know how small the grid would need to be to make the solve time reasonable. And your definition of reasonable would surely differ from mine, depending on how badly you needed the answer.
Javier (view profile)
Hey John, would gridfit allow me to constraint the slope of the fitted surface for a given range of data?
I'm needing to model a surface in the form of z(x,y) from scattered data points while keeping the first derivative of z negative in both directions.
Kim (view profile)
I'm trying to make a surface approximation to describe the required input torque to drive a hydraulic motor operating under certain load conditions. My inputs are: output flow, output pressure and input shaft speed. I have scattered measurement data where the relationship between these variables are to be found. I read in some of the earlier posts that gridfit could be extended to higherdimension fits. Do you have any version where a 3rd dimension can be included that I could try?
Jamal (view profile)
Gridfit is very nice. I have written an update that adds bicubic interpolation and improves the internal scaling (makes surface smoothness completely independent of the x/y grid size). If I send you the code, would you update Gridfit?
Jon Griffiths (view profile)
John D'Errico (view profile)
Alejandro  while I would like to implement that capability one day in the future, that day may be a long time away before I find the time. Until then, all I can offer is to use a finer grid, which will mimic a higher order interpolant.
John D'Errico (view profile)
Dominik  Sorry, but no. You could use this tool to generate a gridded surface, then use interp2, to interpolate, but that is perhaps too much if you only wanted to interpolate a few scattered points.
dominik (view profile)
Hey
i wonder, can this code also be used like an interpolant (e.g. gridded interpolant), i.e. such that the output is not the smooth funtion on a grid, but on certain points queried?
Alejandro Arrizabalaga (view profile)
Great function!
I was just wondering: would it be difficult to extend it for a cubic interpolation? (it was mentioned in the code as a future enhancement).
Roman Voronov (view profile)
Roman Voronov (view profile)
Exactly what I was looking for, the help could start with a simpler example:
[X,Y] = meshgrid(1:5);
Zsmooth = gridfit(X,Y,A,1:5,1:5,'smooth',10)
Zsmooth =
0.9638 1.9517 2.9512 3.9517 4.9638
1.9517 2.8162 3.7338 4.8162 5.9517
2.9512 3.7338 4.5266 5.7338 6.9512
3.9517 4.8162 5.7338 6.8162 7.9517
4.9638 5.9517 6.9512 7.9517 8.9638
John D'Errico (view profile)
Andrew  It sounds like your goal is to reproduce the local shape of your surface over a region that lies inside the bounds of your data. Gridfit can do this, if you are willing to expand that grid so that it fully contains the data. Simply add some nodes that extend to the boundaries of your data in both x and y. The nodes used by gridfit need not be equally spaced, so this should not be a problem, but to retain the essential shape of the surface, you may need to add a few such nodes in each direction. Once the fit is obtained, then you can always extract only that portion of the surface that you wish to retain.
If this is not what you are looking for, then feel free to contact me directly, perhaps with some data so that I can better understand your problem.
Andrew (view profile)
This is an amazing program. I do have one question. I have a scattered grid of latitude(Y) and longitude(X) coordinates corresponding to brightness(Z) values. I've been using GRIDDATA to interpolate to a new grid that is not scattered and within a smaller range of latitude and longitude values, but GRIDDATA smooths the data so much that it becomes useless.
I understand that GRIDFIT is not an interpolant, and from my reading of the help, it does not seem that it was meant to be used in this fashion, but I just want to make sure.
John D'Errico (view profile)
Collin  What you wish is essentially a 3d version of gridfit, that could use time as the third dimension. Sorry, but perhaps one day...
John D'Errico (view profile)
Atul Ingle  Yes, it is possible to extend gridfit to higher dimensions. In practice it tends to turn into a memory & cpu hog before long. 5 or 6 dimensions were the practical extreme when I've done it before, although computers are bigger and faster since those days. 64 bit MATLAB and lots of RAM will help here, but even with that, the curse of dimensionality is huge and arrives quickly. Sadly, the nd version I wrote many years ago is not mine to give away.
Atul Ingle (view profile)
Is it possible to extend this idea to fit 3d volume to scattered data, and in general, to higher dimensions? The only limitation I can see is matrix size.
Collin (view profile)
Excellent bit of code that has made my life a lot easier.
If you have knowledge of the evolution of your surface in time, how would you suggest incorporating knowledge of a previous state to help improve interpolation/extrapolation of the current state?
Thanks
Collin
John D'Errico (view profile)
Laura, I do have a tool (based on the same philosophy as gridfit) to handle fitting over a nonrectangular domain. That set of tools can also test for points inside a general domain, although inpolygon can handle that problem in 2dimensions. I've not posted the code although it is fully functional, because I've never been sufficiently happy with the entire toolset. feel free to send me email if you want to try it though.
Laura (view profile)
I recently (today) submitted a question entitled "availble solutions for 2D interpolation on non monotonic scattered data". Right after I discovered gridfit.
Thanks a lot for such a routine, really solid and with the best fitting I was looking for!.
However I would like to ask one thing that partially has already been addressed. Indeed I am most interested in the good interpolation it performs (although not the objective) than the extrapolation. My x y data points do not define a square domain but more a distorted "trapezoidal" one with curved contours. I would like to get rid of the data out of the original domain.
Ignacio proposed an interesting way (maybe not very clean but working) to know which points fall inside or outside the domain. But I can not apply it to my case since due to these convex borders the Delaunay triangulation produces spurious triangles there that indeed are out of the domain.
Could you suggest me in which way I could identify which of the points of my "square box grid" fall actually inside the trapezoidal domain?
Thanks a lot
John D'Errico (view profile)
Bruno  There are a couple of issues to watch out for here.
First of all, if you wish to find the area that is strictly above zero and below the surface, be careful if the surface goes below zero. (Some individuals miss this point, but it is crucial.) You could simply compute the area by calling trapz(trapz(max(0,Z))). This would be an approximation of course, missing out on the subtle issue of exactly where the surface passes below zero.
Next, depending on how the surface is intended to be interpolated, if you want the EXACT integral of that volume, then be careful. For the first order tensor product surface, i.e., that which interp2 would call 'linear', or what is called 'bilinear' in gridfit, then it suffices to interpolate at the center of each cell in the lattice of points. Sum up those values and you will have an exact volume. (Be careful to scale this result by the width of the cells in each dimension.) As pointed out, trapz (called twice, so once in each dimension) will also be exact here. It turns out that merely computing the average value over each of these cells is also equivalent to that act, then sum up the volumes in each cell of the lattice.
For the triangular faceted surface, trapz is not exact, but the difference may be subtle. Here we end up with a weighted linear combination of the values of each corner of the square cell in the lattice, but the 4 corners are not weighted equally. Send me an email and I can detail how to do this more clearly.
Really though, the easy answer is to use trapz, but you should not ignore these subtleties I've pointed out.
IainIsm (view profile)
John  thank you, 'springs' seems to be doing the job!
Bruno  you can use trapz twice: http://blogs.mathworks.com/loren/2011/06/13/calculatingtheareaunderasurface/
Bruno (view profile)
Excellent work. I have one question... how can I calculate the volume beneath the interpolated surface and above the plane xy? Maybe it's a simple task but i can't figure out how it's done.
John D'Errico (view profile)
IainIsm  You cannot (easily) place bound constraints on this fit. That would result in a very large scale bound constrained fit. While doable in theory with a quadratic programming tool or lsqlin, these problems are huge. It would simply take too long to solve.
Having said that, there is a possible solution that might work. Choose the 'springs' regularizer as an option. That option was designed to avoid extrapolation wherever possible, so if your data never goes negative, then the surface will hopefully not do so either.
John
IainIsm (view profile)
Lovely function  just one question, is it possible to put constraints on the range of values that zgrid can take please? for example, I'm using this to create a surface for a motor rpm/torque/efficiency map (which I am then querying with interp2 to provide an efficiency estimate for a given speed and loading), but every now and then I find that some of the zgrid values are negative. I could set these to zero using find, but I'm hoping to find(!) something a little more elegant!
Chenjie (view profile)
Rijk (view profile)
Jesse (view profile)
Jordan Malof (view profile)
Great stuff! Additionally, and refreshingly, it worked "right out of the box". Thanks John!
bouncing (view profile)
Josh (view profile)
John,
Thanks a lot. Not only did I enjoy an excellent bit of code, but also for your variety of stimulating and entertaining responses. Thank you for breaking up my day with at least 5 minutes of enjoyable rhetoric and exceptional references (Mark Twain / god and little green apples) ... brilliant! Have a few extra stars.
Josh
K E (view profile)
Lovely, thanks
Rodrigo R. Oliveira (view profile)
It is wonderfull. Easy solve my problem!
Gonzalo (view profile)
Peter (view profile)
Oops. In my explanation above P = 1/sqrt(2), not sqrt(2).
Peter (view profile)
Lanis, John, et al.
This is how to get consistent behavior with different grid sizes.

Save a copy of gridfit.m to gridfit2.m.
In gridfit2.m make the following changes at the locations of the commentedout code (near Lines 550600):
% Sorry, John D.
% dy1 = dy(j(:)1)/params.yscale;
% dy2 = dy(j(:))/params.yscale;
dy1 = dy(j(:)1);
dy2 = dy(j(:));
% Sorry, John D.
% Append the regularizer to the interpolation equations,
% scaling the problem first. Use the 1norm for speed.
%NA = norm(A,1);
%NR = norm(Areg,1);
%A = [A;Areg*(smoothparam*NA/NR)];
FidelityEquationCount = size(A, 1);
RegularizerEquationCount = nx * (ny  2) + ny * (nx  2);
NewSmoothParam = smoothparam * sqrt(FidelityEquationCount / RegularizerEquationCount);
A = [A; Areg * NewSmoothParam];

Run the following script to demonstrate:
x = [0, 1.333333333, 2.666666667, 4, 0, 0, 4, 4, 2];
y = [2, 2, 2.1, 1.8114, 0, 4, 0, 4, 0];
z = [0.5, 0.82, 0.4, 0.7, 0.5, 0.7, 0, 0.7, 0];
% Set the smoothness for the original GridFit. This smoothness produces inconsistent results if the grid size changes.
Smoothness = 10;
% Set up the grids with different spacings...
x_coarse = linspace(0, 4, 11);
y_coarse = x_coarse;
x_fine = linspace(0, 4, 41);
y_fine = x_fine;
% We're taking a slice through each surface at the same coordinate.
% They have different indexes but the same coordinate of 1.2.
Index_coarse = find(y_coarse == 1.2);
Index_fine = find(y_fine == 1.2);
% Run GridFit.
z_coarse = gridfit(x, y, z, x_coarse, y_coarse, 'Smoothness', Smoothness);
z_fine = gridfit(x, y, z, x_fine, y_fine, 'Smoothness', Smoothness);
% Try the same thing with the updated GridFit. Note that the smoothness parameter is on a different scale compared to
% the original GridFit. However, all choices of smoothness value should produce consistent results for any grid size.
UpdatedGridFitSmoothness = 1;
% Run the updated version of GridFit, GridFit2.
z_coarse2 = gridfit2(x, y, z, x_coarse, y_coarse, 'interp', 'bilinear', 'Smoothness', UpdatedGridFitSmoothness);
z_fine2 = gridfit2(x, y, z, x_fine, y_fine, 'interp', 'bilinear', 'Smoothness', UpdatedGridFitSmoothness);
% Show the results for GridFit...
% These lines only line up when you use a very large or very small smoothness value.
subplot(1, 2, 1)
hold on
plot(x_coarse, z_coarse(Index_coarse, :), 'color', 'blue');
plot(x_fine, z_fine(Index_fine, :), 'color', 'red')
title('Profiles from the Original GridFit');
axis([1, 4, 0.3, 0.6]);
hold off
% Show the results for GridFit2...
% These lines have the same profile for all smoothness values.
subplot(1, 2, 2)
hold on
plot(x_coarse, z_coarse2(Index_coarse, :), 'color', 'blue');
plot(x_fine, z_fine2(Index_fine, :), 'color', 'red')
title('Profiles from the Updated GridFit');
axis([1, 4, 0.3, 0.6]);
hold off
Peter (view profile)
Lanis, John, et al.
How it works:
GridFit is a function that regularizes a dataset. It constructs a series of linear equations that are the "fidelity equations" (how closely the output matches your input data) and the "smoothness equations" (how close the surface's second derivatives are to zero). It assembles these equations and solves them in a way that minimizes the sum of the squared error (SSE) of the output surface.
The key is to recognize that the SSE is just a sum of squared residuals, the fidelity squared residuals (SSE_F) and the smoothness squared residuals (SSE_Sm).
SSE = SSE_F + SSE_Sm
If we have a certain number of input points, that tells us how many terms are in SSE_F. The number of grid points (more precisely, the number of second derivatives within the grid) tells us how many terms are in SSE_Sm.
Imagine that we are trying to fit this parabola: y = x^2 to our data. The fidelity terms are some finite value (which doesn't matter for this example), such as 100. Because we're using this parabola, the second derivative is always 2. If we have a grid such that there are 20 smoothness terms, then SSE_Sm = 20 * (2^2) = 80.
The ratio SSE_F / SSE_Sm = 100 / 80 = 1.25. This ratio determines how to weigh the smoothness against the goodness of fit.
Now suppose that we suddenly double the number of grid points. Now SSE_Sm = 2 * 20 * (2^2) = 160. That means SSE_F / SSE_Sm = 100 / 160 = 0.625. UhOh! The ratio has changed, which means our fit will favor smoothness more than it did before!
The way to fix this is to reduce the values of the second derivatives by multiplying them by a constant, which we'll call P. In this example, if P = sqrt(2), then the individual terms in SSE_Sm will have half of the value they had before. That's ok because there are twice as many of them. In other words, 0.5 * 2 = 1.
Now using P, SSE_Sm = 2 * 20 * ((P * 2) ^ 2) = 80. This provides the same balance of fidelity vs. smoothness that we originally had.
What I did in GridFit2 is a little more complicated because it also considers the smoothness parameter and the number of fidelity equations. I hope this helps!
John D'Errico (view profile)
Peter  if you can do better, feel free to write your own code. (It is NOT the use of the 1norm that impacts your question.)
Peter (view profile)
I'm not convinced that the smoothness parameter produces consistent results for different grid sizes. I think the problem comes from the use of the 1norm instead of just using the numbers of fidelity equations vs. derivatives.
Here is an example.
x = [0, 1.333333333, 2.666666667, 4, 0, 0, 4, 4, 2];
y = [2, 2, 2.1, 1.8114, 0, 4, 0, 4, 0];
z = [0.5, 0.82, 0.4, 0.7, 0.5, 0.7, 0, 0.7, 0];
Smoothness = 10;
x_coarse = linspace(0, 4, 11);
y_coarse = x_coarse;
Index_coarse = find(y_coarse == 1.2);
x_fine = linspace(0, 4, 41);
y_fine = x_fine;
Index_fine = find(y_fine == 1.2);
z_coarse = gridfit(x, y, z, x_coarse, y_coarse, 'Smoothness', Smoothness);
z_fine = gridfit(x, y, z, x_fine, y_fine, 'Smoothness', Smoothness);
hold on
plot(x_coarse, z_coarse(Index_coarse, :), 'color', 'blue');
plot(x_fine, z_fine(Index_fine, :), 'color', 'red')
hold off
Daniel (view profile)
yes, Hendawi Mohamed I removed points below the edge of my data. My data is roughly bounded by a triangle in XY so I used:
gf=gridfit(....);
Z=griddata(... same limits);
N=find(isnan(Z));
gf(N)=NaN;
to removed the complete estimated results. When I plotted my data, the points on the graph corresponded to area were I had data.
John D'Errico (view profile)
You would normally first decide what size result you wish to see, thus the number of nodes in x and y. Then you would choose some amount of smoothing.
I've set it up so that as well as I could so that changing the grid size for a fixed smoothness will have a minimal impact on the overall look of the surface. For example, try this example:
xy = rand(100,2);
z = rand(100,1);
zgrid50 = gridfit(xy(:,1),xy(:,2),z,50,50);
zgrid200 = gridfit(xy(:,1),xy(:,2),z,200,200);
subplot(1,2,1)
surf(zgrid50)
shading interp
subplot(1,2,2)
surf(zgrid200)
shading interp
The two surfaces will look qualitatively similar, although the second one will look less jagged because it is 4 times finer in each dimension. As expected, this is what you gain by using a finer grid, at some cost in time for the solve. In this sense, smoothness is fully decoupled (as well as I could do that) from the size of the grid itself, for any given set of data.
Given a choice of grid, a smoothness of 1 means that there is an equal amount of importance associated with reducing the residuals compared to the overall smoothness of the surface. Be careful however, in making the smoothness too large or too small, as then you may find numerical artifacts generated from the solution. Floating point arithmetic only goes so far.
By the way, a good rule of thumb is if you cannot get the result you wish with a smoothness somewhere between 0.001 and 1000, then I would look to see if there may be a fundamental problem in the data.
As well, if one chose to fit a 2x2 surface through 1000 data points, you will never get a perfect fit regardless of how small a smoothing parameter you choose, UNLESS the data happens to perfectly fit the implicit model used to interpolate inside a unit square. (For example, tensor product linear interpolation uses an implicit local model of a0 + a1*x + a2*y + a3*x*y within any cell.) There simply are not sufficient degrees of freedom in the model for a very coarse grid if you have many data points.
Ianis (view profile)
Hi! Thank you for this useful function ! I needed to fit scattered points with noise added, and Gridfit handles that nicely ! Now, I would like to test with different values of smoothness and number of xy nodes in order to find the best couple of parameters for my problem. So, I have some questions about the 'smoothness' parameter. You explain in the help that the smoothness parameter is ''normalized'' so a value of '1' may generate reasonable results. I would like to know with respect to which parameter do you normalize the ''smoothness'' value ? If I use either 'gradient ' or ' laplacian' as regularizer, is the smoothness defined in the same way ?
Thanks.
Ianis (view profile)
John D'Errico (view profile)
Gridfit is not a tool to do 2d interpolation. Use interp2 instead.
MAN CHIA (view profile)
How could I use gridfit to interpolate a set of 2D data on a image?
Just let all z equals 0 ?
MAN CHIA (view profile)
R (view profile)
John,
All your submissions I have tried are great quality! Thank you very much. I have a directory for your functions, they could easily be part of Matlab's distribution in my opinion.
Re gridfit, do you still have in mind submitting a 3D/ND version of it?
Rodrigo.
Katie (view profile)
Tomas (view profile)
Great work!
John D'Errico (view profile)
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.
Yan (view profile)
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 (view profile)
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.
Peter (view profile)
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 (view profile)
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.
Mark (view profile)
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?
Mark (view profile)
Darin (view profile)
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?
Christian (view profile)
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.
John D'Errico (view profile)
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 (view profile)
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.
M. C Ertem (view profile)
Thank you!!! Saved the day. (Actually, night.)
Yang (view profile)
Hi John,
I tried some smaller parameters (say 0.1 for smoothness), but then I get wavelike 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!
John D'Errico (view profile)
Smaller values of the smoothing parameter give less smoothing, although zero as a value may result in a singularity.
Yang (view profile)
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!
Qiang (view profile)
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
stefano pasquali (view profile)
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 ?
John D'Errico (view profile)
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 2dimensional 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.
Melissa (view profile)
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 Reinert Nash (view profile)
Excellent tool that much improved my visualization of some photographic exposure data. Thanks, John, great to use one of your tools (again).
John Reinert Nash (view profile)
Januka (view profile)
Great piece of work John. No hassle and easy to use. Thanks!
Ignacio (view profile)
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
jose tapia (view profile)
Thank's for your answer John D'Errico.
John D'Errico (view profile)
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 (view profile)
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 2D table lookup using interpolationextrapolation 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.
Mikkel Lund (view profile)
Very nice.. Just what I have been looking for
John D'Errico (view profile)
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 fortytwo 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 fishingpole. And by the same token any person can see that seven hundred and fortytwo years from now the Lower Mississippi will be only a mile and threequarters 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.
Jason (view profile)
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 (view profile)
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.
Lennert (view profile)
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!
Michael Brown (view profile)
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 PT grid. It was a relief to be able to focus on the science rather than on debugging the scheme to process the data.
Toni (view profile)
Sorry!
Still stumbling about these issues.
Thanks for the right directions!
John D'Errico (view profile)
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 (view profile)
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);
Lars Robert (view profile)
excellent job on meteorological balloon data!
Alexis (view profile)
+1 vote for gridfitn ;)
Ben (view profile)
Thanks John for your code. I enjoy it very much!
xfactor y (view profile)
Adam Chapman (view profile)
the new setting "'smoothness',[xsmooth ysmooth]" is a Godsend. Many Thanks
John D'Errico (view profile)
I like Adam's idea, and this is easily enough put into the code.
Adam Chapman (view profile)
This is brilliant, but a capability to vary the degree of smoothing in x and y exclusively would be very useful
Jeff (view profile)
John D'Errico (view profile)
I've responded via direct email, 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.
Janine Bijsterbosch (view profile)
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 (view profile)
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
John D'Errico (view profile)
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 2dimensions of course.
US Patent # 4992861, 4941039
At least for me, actually reading a patent is about as miserable a pastime as I can imagine.
Matteo Niccoli (view profile)
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.
Jeff Evans (view profile)
Panagiotis (view profile)
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!
John D'Errico (view profile)
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 (view profile)
Is it possible to totally avoid any extrapolation?
Félix Martin (view profile)
Christine A. (view profile)
Chamane (view profile)
excellent
Ulrich (view profile)
Michael Adsetts Edberg Hansen (view profile)
John D'Errico (view profile)
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.
Yaroslav Bulatov (view profile)
What is the advantage of using gridfit instead of griddata followed by smoothing?
Daniel Arevalo (view profile)
Hello John, the function is really nice but I'm wondering if it would be possible to use it with nonrectangular domains where the nodes are just specified with (xn,yn) pairs n=1,..., N
Thanks for any help !
Angel Atanasov (view profile)
Thank you mr. D'Errico, for making my life easier.You piece of code is really good.
Chris Sherwood (view profile)
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.
Karl (view profile)
This is spectacular. I hope the Mathworks is paying you some sort of royalty for your efforts!
Michael Dupin (view profile)
Bryan Raines (view profile)
J Ú (view profile)
Michael Jordan (view profile)
Jveer (view profile)
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)?
Alexis (view profile)
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...
Cristina Velasco (view profile)
Versatile smoothing function.
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.
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
it is very good toolbox ?thanks
Excellent utility!! A default setting did everything I needed with minimal error.
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...
1. In strict accordance with the help, xi and yi WERE xnodes and ynodes 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
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.
What is wrong with calling GRIDFIT? Needs improvement.
Z=rand(100);
a=1:100;
[X,Y]=meshgrid(a,a);
good=(X50).^2+(Y50).^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
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
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
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?
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.
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!
Great tool and insightful docs! Gridfit succeeds where griddata fails; that is, with noisy & illspaced data such as exist in the real world. Thanks!
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
This is very very good work!
This is simply great, with all we can expect for a replacement to griddata (which fails too often).
Very very good John! Take good care of Amy!
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 goodlooking 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!
Excellent!!!
Excellent! This is exactly what I was searching for. Thank you very much for this great contribution.
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.
I'll invite Håvard Torpe to contact me via email 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.
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?
Marvelous!
Should be included into the default distribution of matlab! perhaps somehow integrated into griddata?
Does a great job!
Excellent Work!
This is an excellent program. It makes it easy to fit surfaces to 3D data.
Great ! Thank you!
Thank you...
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
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
This is exactly what I am looking for... Thanks you!!! Save me a lot of time to write one...
simply said: just another extremely useful member of the d'erriconiandraconian family of musthaves if you're professionally serious about nddata investigation / reconstruction / interpolation / visualization, which smoothly unites with its nowfamous siblings (consolidator, inpaintnans) into a very nice toolbox for those everyday cheatingwithdata endeavors...
the code is very clean and the profiler reveals no unnecessary hotspot...
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
Neat picture!