File Exchange

image thumbnail

Surface Fitting using gridfit

version (2.56 MB) by John D'Errico
Model 2-d surfaces from scattered data


Updated 04 Mar 2016

View License

Editor's Note: Popular File 2008

This file was selected as MATLAB Central Pick of the Week

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

Cite As

John D'Errico (2021). Surface Fitting using gridfit (, MATLAB Central File Exchange. Retrieved .

Comments and Ratings (278)

steve didienne

combustion turbulence

patel mukesh

Great,, I was so stressed looking at my scattered data and it was not possible to sort them grid them and make it usable to feed in for interpolation methods

Sanjay Manohar

Amazing. Spent hours looking for a solution.

Perveen Sethi

This product seems to perform well relative to your CFC. Thanks again.

Aritra Ghosal

Onur Kurum

qingbin meng


Thanks A lot ! Robust Code with all options you could dream off. Just used it on a set of 50GB of .xyz data splitted into 50 files depending on the discretization of the output you want it runs as fast or even faster than any Geographical Information Software (GIS) on the market. I added a time estimate to completion to help the user visualizing the time needed for very large data set. It would be great to have a code that estimate the best Tile Size for a give xyz. Right now I am fixing the tile to 10% of the biggest size of element x or y with a maximum to 1000.

To be added in line 1019
time_estimate_x = (cputime-tx)* (nx-xtind(end))/tilesize;
disp(['Time Estimate :' num2str(time_estimate_x/3600) 'Hr /' num2str(time_estimate_x/60) 'Min']);

Rashid Yagoub


Massimo Del Guasta

Excellent function! it solved many troubles I had. Thanks a lot!

Ted Shultz

This is more than just a useful function. The code is extremely well documented, and includes a separate documentation file explaining the idea behind the function! I love that this function is both useful for my application, but also that I can learn from the author by reading the notes and looking at the code!




Derya Akkaynak


Jason Nicholson

Masoom Kumar, when I get a chance this summer, I will upload more about how to do this with all the various types of constraints.

Masoom Kumar

Is there any possibility to apply bounds on the fitting ? In other words can we restrict the fit between some upper and lower bounds ?

Jaco Verster

Alessandro Juri

ma os

Omar Marello

Guangpeng Yan

uma shankar

uma shankar

Mahdi S. Hosseini

A relevant toolbox for numerical differentiation called MaxPol has been recently published and made available here

MaxPol provides a framework to design variety of numerical differentiation kernels with properties like:
(1) Cutoff (lowpass) design with no side-lob artifacts (for noise-robust 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


Howerver , I don't understand the difference between this one and "scatteredinterpolant" in the function hub of maltab

Gert Kruger

Jason Nicholson

If you like this function, you should try regularizeNd. It is nD version of gridfit with some bug fixes and improvements.

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

Arun Krishnadas

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

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 sfit-object 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


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 gaussian-like 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

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

Excellent work!


Somayeh Hesabi

pengfa zhou

hi,i'm a learner of matlab, anyone can tell me where i should paste this file for use?thanks million time!

j b

talyn dong



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

rajyalakshmi k

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

it is not working for non-monotone increasing vectors

Behailu Shikur


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!


Lyra Maxwell

Julia Howe

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?



Zhang Lidong

Very Good

Suzie Oman

Does this create a formula of the fit?

Sathish Sanjeevi

yousaf obaid


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


I use it often, thanks for your job



Rob Campbell

daniel perez

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

daniel perez

Amazing work!


Wrichik Basu


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?


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.



Can someone help with this error?
Undefined function 'sparse' for input arguments of type 'single'. Thanks.



Great submission as always John. Thanks.

Already asked below, but can you add bi-cubic interpolation to this submission?

Please :-)





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?


Okay, thank you.

John D'Errico

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.


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

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 n-dimensional version of this long ago, but that code belongs to my past employer.

Jarrod Collins

This code is wonderful.I was wondering, is there a C++ or VTK implementation to this?

John D'Errico

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

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.


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

I have used this program for years and it works well! Has anyone tried converting gridfit to C using Matlab Coder ?

John D'Errico

Gridfit is not an interpolation tool. If you want that, use tools like griddata, triscatteredInterp, scatteredInterpolant, etc.

Jia Le Ngai

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

John Chen

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


Thanks so much. I am very appreciative.


John D'Errico

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.


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?)


John D'Errico


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.


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 1E-10, 1E-5, 1E-1, 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.


John D'Errico

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


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.



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 off-the-wall needs have come up, it's been a great base to start from. Clear and well structured code.

John D'Errico

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

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 wishy-washy warning, that may or may not generate a warning when it encounters something surprising.

Chad Greene

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)])
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)


This is a fantastic function and it runs perfectly with my 2D-data. Do you also provide a gridfit3d version? That whould save a lot of time in my post-processing work. Thanks in advance.


How to cited your code in research paper?

John D'Errico

Evgheny - If you seriously need to fit such a surface in spherical coordinates, I might be able to help you out with a separate, un-posted toolbox that is capable of fitting a surface in spherical coordinates. You would need to contact me directly.


Thanks, function does do its job great!

Does anyone know alternative of this function for closed surface (sphere-like).
If I try this function with spherical coordinates, it works fine, but there are problems with the seam (line of joint)

John D'Errico

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.


is there any implementation of this in java? that get some points location and value and give a grdi?


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.


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!


I found that it is very straightforward to fit a nice looking surface to some 55 data-points 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.


John D'Errico

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 multi-valued, 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

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

Still it cannot surface a cluster of points that looks like a ball. The result the interpolation is always a surface of single-value function.

Highman Field

Angelo Tafuni

John D'Errico

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.


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.


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 higher-dimension fits. Do you have any version where a 3rd dimension can be included that I could try?


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

John D'Errico

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

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.


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?

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

Roman Voronov

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

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.


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

Collin - What you wish is essentially a 3-d version of gridfit, that could use time as the third dimension. Sorry, but perhaps one day...

John D'Errico

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 n-d version I wrote many years ago is not mine to give away.

Atul Ingle

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.


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?



John D'Errico

Laura, I do have a tool (based on the same philosophy as gridfit) to handle fitting over a non-rectangular domain. That set of tools can also test for points inside a general domain, although inpolygon can handle that problem in 2-dimensions. 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 e-mail if you want to try it though.


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

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


John - thank you, 'springs' seems to be doing the job!

Bruno - you can use trapz twice:


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

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.



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!




Jordan Malof

Great stuff! Additionally, and refreshingly, it worked "right out of the box". Thanks 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.



Lovely, thanks

Rodrigo R. Oliveira

It is wonderfull. Easy solve my problem!



Oops. In my explanation above P = 1/sqrt(2), not sqrt(2).


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 commented-out code (near Lines 550-600):

% 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 1-norm 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


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

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. Uh-Oh! 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

Peter - if you can do better, feel free to write your own code. (It is NOT the use of the 1-norm that impacts your question.)


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 1-norm 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


yes, Hendawi Mohamed I removed points below the edge of my data. My data is roughly bounded by a triangle in X-Y so I used:

Z=griddata(... same limits);


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

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);
shading interp
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.


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 x-y- 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 ?


John D'Errico

Gridfit is not a tool to do 2-d interpolation. Use interp2 instead.


How could I use gridfit to interpolate a set of 2D data on a image?
Just let all z equals 0 ?




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?




Great work!

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.


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

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.


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);

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

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.


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?



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?


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

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?


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

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


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!

John D'Errico

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


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


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);


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 ?

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.

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


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

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


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


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

Thank's for your answer John D'Errico.

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

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.

Mikkel Lund

Very nice.. Just what I have been looking for

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.


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

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.


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

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.


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

John D'Errico


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
k = isnan(x) | isnan(y) | isnan(z);
if any(k)

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?


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);
x = 1:512;
y = 1:512;
zgd = griddata(x,y,z,xg,yg);

Lars Robert

excellent job on meteorological balloon data!


+1 vote for gridfitn ;)


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

xfactor y

Adam Chapman

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

John D'Errico

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

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

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

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

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.

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.

Jeff Evans


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

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.


Is it possible to totally avoid any extrapolation?

Félix Martin

Christine A.




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.

Yaroslav Bulatov

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

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 !

Angel Atanasov

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

Chris Sherwood

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


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

Michael Dupin

Bryan Raines


Michael Jordan


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)?


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

Pauline Wong

Versatile smoothing function.

Cedric Testaz

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.

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 . Compare with your own method. My own is in fact interpolation, but could be modified for approximation without much problems.
Andrey, Kamchatka

gai cs

it is very good toolbox ?thanks

Anmol Agrawal

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

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

Alex Gardner


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.

ans =

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.


What is wrong with calling GRIDFIT? Needs improvement.

??? Error using ==> gridfit at 404
xnodes and ynodes must be monotone increasing

Tim Marvin

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);

Thank you,

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

Henrique Mendes

It is not working! Matlab says:
"??? Input argument "x" is undefined.
Error in ==> gridfit at 325
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?

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.

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!

Juan Perez

gabriele flauti

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!

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

Zhijun Wang

This is very very good work!

E Farhi

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

Carlos Adrián Vargas Aguilera

Very very good John! Take good care of Amy!

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!

Phani Ivatury


Wolfgang Schwanghart

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

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.

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.

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?

Andy Lulham

thank you


Aslak Grinsted

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

Wen Soong

Does a great job!


Excellent Work!

Kevin Denis

David Barker

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

Julio Oliveira

Great ! Thank you!


Thank 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

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

Chaowen Liang

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

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

Michael Robbins

Neat picture!

MATLAB Release Compatibility
Created with R2016a
Compatible with any release
Platform Compatibility
Windows macOS Linux

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!