(This will get long, unless I just run out of steam.)
First, a few comments about non-scattered data. (See the answer that references littleCMS)
There are two types of color interpolation that is common. Some years ago, trilinear interpolation (tensor product linear interpolation) was a common approach for color table interpolation. Trilinear interpolation can indeed be implemented as a sequential set of one dimensional interpolations, first on one axis, then then along a second axis, etc.
Many years ago we all realized that trilinear interpolation introduces artifacts in color imaging, when applied to certain types of transformations. The problems are seen in neutrals. A solution is to move to a simplicial interpolant, in 3-d, by dissecting a cube into 6 tetrahedra. In n dimensions, the unit cube will be dissected into factorial(n) simplexes. There are other dissections of a cube, but this particular style assures that the main diagonal is always a shared edge for all the simplexes. This in turn restores good behavior for neutrals, when applied to certain color lookup tables.
Now let me get into the question of true scattered data interpolation.
Others have mentioned a variety of schemes. Kriging, multiquadrics, distance based methods are a few. (When I did some work in the past with these schemes, I actually preferred the inverse multiquadric methods.) All of these are really just variations of radial basis function methods, a common scheme. RBF methods have their good and bad points. They will usually generate a smooth interpolant, this of course depends on the specific basis function chosen, as well as whether you choose to limit the support. RBF methods also allow you to extrapolate, at least as far out as the support of the radial basis elements will extend. If the basis elements are allowed to be infinite in extent, then no explicit constraint on extrapolation will apply. (Extrapolation in general is a bad thing to do.) One problem with RBF methods is they require the solution of large systems of linear equations, and those equation systems are often dense matrices. This means the size of the problem, in terms of the number of data points you can handle tends to be limited by the linear algebra. If instead, you limit the support by truncating the basis elements, then the matrices can become sparse. This will improve the linear algebra if you use a sparse matrix package for the solution. At the same time, the support distance becomes a nonlinear parameter that you must control. As well, methods like multiquadrics and inverse multiquadric methods may have a secondary nonlinear parameter that controls the shape of the basis elements. Kriging has similar issues, and I'd lump all of these methods together.
Do to these issues, all of these methods that I've classed as RBF variants are often limited in the number of points they will comfortably handle. Depending upon how you deal with things and the amount of memory available, that limit might often be on the order of a few thousand points.
One other problem with the general class of RBF methods is what I'll call intrapolation. This is a neologism I created many years ago to describe interpolation across a relatively large hole in the data. In fact, there may often be problems even when interpolating across smaller holes in the data. These methods, because they are smooth to some extent, may introduce unwanted extrema (large peaks or valleys) into the interpolated surface. This is a common problem with even 1-d interpolants, often seen as ringing artifacts with cubic spline or polynomial interpolants, and certainly seen with Fourier series interpolants. The problem in higher dimensions is to even recognize that it has indeed happened, since plotting surfaces in more than three dimensions tends to be difficult.
If you have more points than that limit, or if these ringing artifacts are unacceptable, then other methods are often a better choice. If you are willing to use a linear interpolant, then the simplest solution in higher dimensions is to start with a tessellation of the data. Thus in 3 dimensions, tessellate the data (typically a delaunay tessellation) into tetrahedra. This is fairly efficient to do, and there are many tools to be found for this purpose. It is then a simple problem to interpolate any individual point. Merely identify which simplex the point lies in, compute barycentric coordinates as interpolation weights within the simplex, and form the corresponding linear combination of the function values at each vertex of the found simplex. This is all extremely fast and efficient.
A downside of these tessellation based methods is they typically restrict you to the convex hull of the data points, and as bad, if your data happens to lie in a non-convex domain, then the interpolant may well do strange things in some regions of your domain. Another problem with the scheme I mentioned above, is the interpolant will only be piecewise linear, but once you move into higher dimensions things get nasty fast. Other methods are to be found for smooth interpolation based on a tessellation, but they will take more effort and are therefore far less common.
The basic tradeoffs should be obvious here. If you need a smooth interpolant and only have a few points, then RBF methods are often chosen. They are simple, easy to use, etc. The actual method chosen is often just a matter of convenience, or even habit. I'f I've used one tool before and was happy, I'll probably be happy with it again. Since the question was which method is "best for practical use", I'll point out that best is a very subjective word when applied out of context. What are your goals in an interpolation problem? What skill set do you own? What set of tools do you know how to use? What environment will you work in? All of these factors will influence your choice of the best method.
If you have many data points, and speed is of the essence, but ultimate smoothness is not that important, then you will generally look for a simplicial interpolant. Of course, if you have sufficient points, then the piecewise linear nature of the beast is of less importance. The piecewise linear interpolant here has the great virtue in some cases that it can never generate extrema in your surface that did not exist in the data. For some problems, color characterization for example, this is of paramount importance.
One other issue is with noise. While the presence of noise is often a signal that smoothing of some sort is necessary, not all such surfaces have smoothing applied. Any smoothing operator will sometimes smooth out important features of the data too. This happens because we can think of a smoothing operator as a low pass filter. High frequency behavior is often noise, but it may also be just a sharp toe or shoulder in my surface that I cannot afford to lose. If this is a problem, then you may be willing to use an interpolant even in the presence of sometimes significant noise. In that event, I'll suggest that the simplest, lowest order interpolant is best. A smooth, more global interpolant will also tend to amplify any noise in the data, so if you look for the lowest variance interpolant in the presence of noise, it will generally be a linear interpolant.
Of course, there are many varieties of thin plate splines, interpolatory or not. Once you go beyond one dimension, your options also expand, at least if you are willing to do the work.
I'll end here before it turns into a book.