function [zgrid,xgrid,ygrid] = gridfit(x,y,z,xnodes,ynodes,varargin)% gridfit: estimates a surface on a 2d grid, based on scattered data%          Replicates are allowed. All methods extrapolate to the grid%          boundaries. Gridfit uses a modified ridge estimator to%          generate the surface, where the bias is toward smoothness.%%          Gridfit is not an interpolant. Its goal is a smooth surface%          that approximates your data, but allows you to control the%          amount of smoothing.%% usage #1: zgrid = gridfit(x,y,z,xnodes,ynodes);% usage #2: [zgrid,xgrid,ygrid] = gridfit(x,y,z,xnodes,ynodes);% usage #3: zgrid = gridfit(x,y,z,xnodes,ynodes,prop,val,prop,val,...);%% Arguments: (input)%  x,y,z - vectors of equal lengths, containing arbitrary scattered data%          The only constraint on x and y is they cannot ALL fall on a%          single line in the x-y plane. Replicate points will be treated%          in a least squares sense.%%          ANY points containing a NaN are ignored in the estimation%%  xnodes - vector defining the nodes in the grid in the independent%          variable (x). xnodes need not be equally spaced. xnodes%          must completely span the data. If they do not, then the%          'extend' property is applied, adjusting the first and last%          nodes to be extended as necessary. See below for a complete%          description of the 'extend' property.%%          If xnodes is a scalar integer, then it specifies the number%          of equally spaced nodes between the min and max of the data.%%  ynodes - vector defining the nodes in the grid in the independent%          variable (y). ynodes need not be equally spaced.%%          If ynodes is a scalar integer, then it specifies the number%          of equally spaced nodes between the min and max of the data.%%          Also see the extend property.%%  Additional arguments follow in the form of property/value pairs.%  Valid properties are:%    'smoothness', 'interp', 'regularizer', 'solver', 'maxiter'%    'extend', 'tilesize', 'overlap'%%  Any UNAMBIGUOUS shortening (even down to a single letter) is%  valid for property names. All properties have default values,%  chosen (I hope) to give a reasonable result out of the box.%%   'smoothness' - scalar - determines the eventual smoothness of the%          estimated surface. A larger value here means the surface%          will be smoother. Smoothness must be a non-negative real%          number.%%          Note: the problem is normalized in advance so that a%          smoothness of 1 MAY generate reasonable results. If you%          find the result is too smooth, then use a smaller value%          for this parameter. Likewise, bumpy surfaces suggest use%          of a larger value. (Sometimes, use of an iterative solver%          with too small a limit on the maximum number of iterations%          will result in non-convergence.)%%          DEFAULT: 1%%%   'interp' - character, denotes the interpolation scheme used%          to interpolate the data.%%          DEFAULT: 'triangle'%%          'bilinear' - use bilinear interpolation within the grid%                     (also known as tensor product linear interpolation)%%          'triangle' - split each cell in the grid into a triangle,%                     then linear interpolation inside each triangle%%          'nearest' - nearest neighbor interpolation. This will%                     rarely be a good choice, but I included it%                     as an option for completeness.%%%   'regularizer' - character flag, denotes the regularization%          paradignm to be used. There are currently three options.%%          DEFAULT: 'gradient'%%          'diffusion' or 'laplacian' - uses a finite difference%              approximation to the Laplacian operator (i.e, del^2).%%              We can think of the surface as a plate, wherein the%              bending rigidity of the plate is specified by the user%              as a number relative to the importance of fidelity to%              the data. A stiffer plate will result in a smoother%              surface overall, but fit the data less well. I've%              modeled a simple plate using the Laplacian, del^2. (A%              projected enhancement is to do a better job with the%              plate equations.)%%              We can also view the regularizer as a diffusion problem,%              where the relative thermal conductivity is supplied.%              Here interpolation is seen as a problem of finding the%              steady temperature profile in an object, given a set of%              points held at a fixed temperature. Extrapolation will%              be linear. Both paradigms are appropriate for a Laplacian%              regularizer.%%          'gradient' - attempts to ensure the gradient is as smooth%              as possible everywhere. Its subtly different from the%              'diffusion' option, in that here the directional%              derivatives are biased to be smooth across cell%              boundaries in the grid.%%              The gradient option uncouples the terms in the Laplacian.%              Think of it as two coupled PDEs instead of one PDE. Why%              are they different at all? The terms in the Laplacian%              can balance each other.%%          'springs' - uses a spring model connecting nodes to each%              other, as well as connecting data points to the nodes%              in the grid. This choice will cause any extrapolation%              to be as constant as possible.%%              Here the smoothing parameter is the relative stiffness%              of the springs connecting the nodes to each other compared%              to the stiffness of a spting connecting the lattice to%              each data point. Since all springs have a rest length%              (length at which the spring has zero potential energy)%              of zero, any extrapolation will be minimized.%%          Note: I don't terribly like the 'springs' strategy.%          It tends to drag the surface towards the mean of all%          the data. Its been left in only because the paradigm%          interests me.%%%   'solver' - character flag - denotes the solver used for the%          resulting linear system. Different solvers will have%          different solution times depending upon the specific%          problem to be solved. Up to a certain size grid, the%          direct \ solver will often be speedy, until memory%          swaps causes problems.%%          What solver should you use? Problems with a significant%          amount of extrapolation should avoid lsqr. \ may be%          best numerically for small smoothnesss parameters and%          high extents of extrapolation.%%          Large numbers of points will slow down the direct%          \, but when applied to the normal equations, \ can be%          quite fast. Since the equations generated by these%          methods will tend to be well conditioned, the normal%          equations are not a bad choice of method to use. Beware%          when a small smoothing parameter is used, since this will%          make the equations less well conditioned.%%          DEFAULT: 'normal'%%          '\' - uses matlab's backslash operator to solve the sparse%                     system. 'backslash' is an alternate name.%%          'symmlq' - uses matlab's iterative symmlq solver%%          'lsqr' - uses matlab's iterative lsqr solver%%          'normal' - uses \ to solve the normal equations.%%%   'maxiter' - only applies to iterative solvers - defines the%          maximum number of iterations for an iterative solver%%          DEFAULT: min(10000,length(xnodes)*length(ynodes))%%%   'extend' - character flag - controls whether the first and last%          nodes in each dimension are allowed to be adjusted to%          bound the data, and whether the user will be warned if%          this was deemed necessary to happen.%%          DEFAULT: 'warning'%%          'warning' - Adjust the first and/or last node in%                     x or y if the nodes do not FULLY contain%                     the data. Issue a warning message to this%                     effect, telling the amount of adjustment%                     applied.%%          'never'  - Issue an error message when the nodes do%                     not absolutely contain the data.%%          'always' - automatically adjust the first and last%                     nodes in each dimension if necessary.%                     No warning is given when this option is set.%%%   'tilesize' - grids which are simply too large to solve for%          in one single estimation step can be built as a set%          of tiles. For example, a 1000x1000 grid will require%          the estimation of 1e6 unknowns. This is likely to%          require more memory (and time) than you have available.%          But if your data is dense enough, then you can model%          it locally using smaller tiles of the grid.%%          My recommendation for a reasonable tilesize is%          roughly 100 to 200. Tiles of this size take only%          a few seconds to solve normally, so the entire grid%          can be modeled in a finite amount of time. The minimum%          tilesize can never be less than 3, although even this%          size tile is so small as to be ridiculous.%%          If your data is so sparse than some tiles contain%          insufficient data to model, then those tiles will%          be left as NaNs.%%          DEFAULT: inf%%%   'overlap' - Tiles in a grid have some overlap, so they%          can minimize any problems along the edge of a tile.%          In this overlapped region, the grid is built using a%          bi-linear combination of the overlapping tiles.%%          The overlap is specified as a fraction of the tile%          size, so an overlap of 0.20 means there will be a 20%%          overlap of successive tiles. I do allow a zero overlap,%          but it must be no more than 1/2.%%          0 <= overlap <= 0.5%%          Overlap is ignored if the tilesize is greater than the%          number of nodes in both directions.%%          DEFAULT: 0.20%%%   'autoscale' - Some data may have widely different scales on%          the respective x and y axes. If this happens, then%          the regularization may experience difficulties. %          %          autoscale = 'on' will cause gridfit to scale the x%          and y node intervals to a unit length. This should%          improve the regularization procedure. The scaling is%          purely internal. %%          autoscale = 'off' will disable automatic scaling%%          DEFAULT: 'on'%%% Arguments: (output)%  zgrid   - (nx,ny) array containing the fitted surface%%  xgrid, ygrid - as returned by meshgrid(xnodes,ynodes)%%% Speed considerations:%  Remember that gridfit must solve a LARGE system of linear%  equations. There will be as many unknowns as the total%  number of nodes in the final lattice. While these equations%  may be sparse, solving a system of 10000 equations may take%  a second or so. Very large problems may benefit from the%  iterative solvers or from tiling.%%% Example usage:%%  x = rand(100,1);%  y = rand(100,1);%  z = exp(x+2*y);%  xnodes = 0:.1:1;%  ynodes = 0:.1:1;%%  g = gridfit(x,y,z,xnodes,ynodes);%% Note: this is equivalent to the following call:%%  g = gridfit(x,y,z,xnodes,ynodes, ...%              'smooth',1, ...%              'interp','triangle', ...%              'solver','normal', ...%              'regularizer','gradient', ...%              'extend','warning', ...%              'tilesize',inf);%%% Author: John D'Errico% e-mail address: woodchips@rochester.rr.com% Release: 2.0% Release date: 5/23/06% set defaultsparams.smoothness = 1;params.interp = 'triangle';params.regularizer = 'gradient';params.solver = 'normal';params.maxiter = [];params.extend = 'warning';params.tilesize = inf;params.overlap = 0.20;params.mask = []; params.autoscale = 'on';params.xscale = 1;params.yscale = 1;% was the params struct supplied?if ~isempty(varargin)  if isstruct(varargin{1})    % params is only supplied if its a call from tiled_gridfit    params = varargin{1};    if length(varargin)>1      % check for any overrides      params = parse_pv_pairs(params,varargin{2:end});    end  else    % check for any overrides of the defaults    params = parse_pv_pairs(params,varargin);  endend% check the parameters for acceptabilityparams = check_params(params);% ensure all of x,y,z,xnodes,ynodes are column vectors,% also drop any NaN datax=x(:);y=y(:);z=z(:);k = isnan(x) | isnan(y) | isnan(z);if any(k)  x(k)=[];  y(k)=[];  z(k)=[];endxmin = min(x);xmax = max(x);ymin = min(y);ymax = max(y);% did they supply a scalar for the nodes?