Home > matlabmk > griddata2.m

# griddata2

## PURPOSE

GRIDDATA Data gridding and surface fitting.

## SYNOPSIS

function [xi,yi,zi] = griddata2(x,y,z,xi,yi,method,options)

## DESCRIPTION

```GRIDDATA Data gridding and surface fitting.

GRIDDATA is not recommended. Use TriScatteredInterp instead.

ZI = GRIDDATA(X,Y,Z,XI,YI) fits a surface of the form Z = F(X,Y) to the
data in the (usually) nonuniformly-spaced vectors (X,Y,Z). GRIDDATA
interpolates this surface at the points specified by (XI,YI) to produce
ZI.  The surface always goes through the data points. XI and YI are
usually a uniform grid (as produced by MESHGRID) and is where GRIDDATA
gets its name.

XI can be a row vector, in which case it specifies a matrix with
constant columns. Similarly, YI can be a column vector and it specifies
a matrix with constant rows.

[XI,YI,ZI] = GRIDDATA(X,Y,Z,XI,YI) also returns the XI and YI formed
this way (the results of [XI,YI] = MESHGRID(XI,YI)).

[...] = GRIDDATA(X,Y,Z,XI,YI,METHOD) where METHOD is one of
'linear'    - Triangle-based linear interpolation (default)
'cubic'     - Triangle-based cubic interpolation
'nearest'   - Nearest neighbor interpolation
'v4'        - MATLAB 4 griddata method
defines the type of surface fit to the data. The 'cubic' and 'v4'
methods produce smooth surfaces while 'linear' and 'nearest' have
discontinuities in the first and zero-th derivative respectively.  All
the methods except 'v4' are based on a Delaunay triangulation of the
data.
If METHOD is [], then the default 'linear' method will be used.

[...] = GRIDDATA(X,Y,Z,XI,YI,METHOD,OPTIONS) specifies a cell array of
strings OPTIONS that were previously used by Qhull. Qhull-specific
OPTIONS are no longer required and are currently ignored. Support for
these options will be removed in a future release.

Example:
x = rand(100,1)*4-2; y = rand(100,1)*4-2; z = x.*exp(-x.^2-y.^2);
ti = -2:.25:2;
[xi,yi] = meshgrid(ti,ti);
zi = griddata(x,y,z,xi,yi);
mesh(xi,yi,zi), hold on, plot3(x,y,z,'o'), hold off

See also TriScatteredInterp, DelaunayTri, GRIDDATAN, DELAUNAY,
INTERP2, MESHGRID, DELAUNAYN.

Note, this is the same as the Matlab function griddata.m save that
"&&~strcmp(method,'invdist')" has been added to the line "if ( nargin < 6
|| isempty(method) ) &&~strcmp(method,'invdist'),  method = 'linear';
end".  This is necessary because in Matlab V12 they changed griddata in
such a way that it is no longer compatible with EEGLAB's topoplot
function (or my topoplotMK.m version of topoplot.m).```

## CROSS-REFERENCE INFORMATION

This function calls:
This function is called by:
• topoplotMK topoplotMK() - plot a topographic map of a scalp data field in a 2-D circular view

## SOURCE CODE

```0001 function [xi,yi,zi] = griddata2(x,y,z,xi,yi,method,options)
0002 %GRIDDATA Data gridding and surface fitting.
0003 %
0004 %   GRIDDATA is not recommended. Use TriScatteredInterp instead.
0005 %
0006 %   ZI = GRIDDATA(X,Y,Z,XI,YI) fits a surface of the form Z = F(X,Y) to the
0007 %   data in the (usually) nonuniformly-spaced vectors (X,Y,Z). GRIDDATA
0008 %   interpolates this surface at the points specified by (XI,YI) to produce
0009 %   ZI.  The surface always goes through the data points. XI and YI are
0010 %   usually a uniform grid (as produced by MESHGRID) and is where GRIDDATA
0011 %   gets its name.
0012 %
0013 %   XI can be a row vector, in which case it specifies a matrix with
0014 %   constant columns. Similarly, YI can be a column vector and it specifies
0015 %   a matrix with constant rows.
0016 %
0017 %   [XI,YI,ZI] = GRIDDATA(X,Y,Z,XI,YI) also returns the XI and YI formed
0018 %   this way (the results of [XI,YI] = MESHGRID(XI,YI)).
0019 %
0020 %   [...] = GRIDDATA(X,Y,Z,XI,YI,METHOD) where METHOD is one of
0021 %       'linear'    - Triangle-based linear interpolation (default)
0022 %       'cubic'     - Triangle-based cubic interpolation
0023 %       'nearest'   - Nearest neighbor interpolation
0024 %       'v4'        - MATLAB 4 griddata method
0025 %   defines the type of surface fit to the data. The 'cubic' and 'v4'
0026 %   methods produce smooth surfaces while 'linear' and 'nearest' have
0027 %   discontinuities in the first and zero-th derivative respectively.  All
0028 %   the methods except 'v4' are based on a Delaunay triangulation of the
0029 %   data.
0030 %   If METHOD is [], then the default 'linear' method will be used.
0031 %
0032 %   [...] = GRIDDATA(X,Y,Z,XI,YI,METHOD,OPTIONS) specifies a cell array of
0033 %   strings OPTIONS that were previously used by Qhull. Qhull-specific
0034 %   OPTIONS are no longer required and are currently ignored. Support for
0035 %   these options will be removed in a future release.
0036 %
0037 %   Example:
0038 %      x = rand(100,1)*4-2; y = rand(100,1)*4-2; z = x.*exp(-x.^2-y.^2);
0039 %      ti = -2:.25:2;
0040 %      [xi,yi] = meshgrid(ti,ti);
0041 %      zi = griddata(x,y,z,xi,yi);
0042 %      mesh(xi,yi,zi), hold on, plot3(x,y,z,'o'), hold off
0043 %
0044 %   See also TriScatteredInterp, DelaunayTri, GRIDDATAN, DELAUNAY,
0045 %   INTERP2, MESHGRID, DELAUNAYN.
0046 %
0047 % Note, this is the same as the Matlab function griddata.m save that
0048 % "&&~strcmp(method,'invdist')" has been added to the line "if ( nargin < 6
0049 % || isempty(method) ) &&~strcmp(method,'invdist'),  method = 'linear';
0050 % end".  This is necessary because in Matlab V12 they changed griddata in
0051 % such a way that it is no longer compatible with EEGLAB's topoplot
0052 % function (or my topoplotMK.m version of topoplot.m).
0053
0054 %   Copyright 1984-2009 The MathWorks, Inc.
0055 %   \$Revision: 5.33.4.11 \$  \$Date: 2009/04/21 03:25:31 \$
0056
0057 error(nargchk(5,7,nargin,'struct'));
0058
0059 [msg,x,y,z,xi,yi] = xyzchk(x,y,z,xi,yi);
0060 if ~isempty(msg), error(msg); end
0061 if ndims(x) > 2 || ndims(y) > 2 || ndims(xi) > 2 || ndims(yi) > 2
0062     error('MATLAB:griddata:HigherDimArray',...
0063           'X,Y and XI,YI cannot be arrays of dimension greater than two.');
0064 end
0065
0066 if ( issparse(x) || issparse(y) || issparse(z) || issparse(xi) || issparse(yi) )
0067     error('MATLAB:griddata:InvalidDataSparse',...
0068         'Input data cannot be sparse.');
0069 end
0070
0071 if ( ~isreal(x) || ~isreal(y) || ~isreal(xi) || ~isreal(yi) )
0072     error('MATLAB:griddata:InvalidDataComplex',...
0073         'Input data cannot be complex.');
0074 end
0075
0076 if ( nargin < 6 || isempty(method) ) &&~strcmp(method,'invdist'),  method = 'linear'; end
0077 if ~ischar(method),
0078   error('MATLAB:griddata:InvalidMethod',...
0079         'METHOD must be one of ''linear'',''cubic'',''nearest'', or ''v4''.');
0080 end
0081
0082 if nargin == 7
0083     if ~iscellstr(options)
0084         error('MATLAB:OptsNotStringCell',...
0085               'OPTIONS should be cell array of strings.');
0086     end
0087     opt = options;
0088 else
0089     opt = [];
0090 end
0091
0092 if numel(x) < 3 || numel(y) < 3
0093   error('MATLAB:griddata:NotEnoughSamplePts',...
0094         'Not enough unique sample points specified.');
0095 end
0096
0097
0098 % Sort x and y so duplicate points can be averaged before passing to delaunay
0099
0100 %Need x,y and z to be column vectors
0101 sz = numel(x);
0102 x = reshape(x,sz,1);
0103 y = reshape(y,sz,1);
0104 z = reshape(z,sz,1);
0105 sxyz = sortrows([x y z],[2 1]);
0106 x = sxyz(:,1);
0107 y = sxyz(:,2);
0108 z = sxyz(:,3);
0109 myepsx = eps(0.5 * (max(x) - min(x)))^(1/3);
0110 myepsy = eps(0.5 * (max(y) - min(y)))^(1/3);
0111 ind = [0; ((abs(diff(y)) < myepsy) & (abs(diff(x)) < myepsx)); 0];
0112
0113 if sum(ind) > 0
0114   warning('MATLAB:griddata:DuplicateDataPoints',['Duplicate x-y data points ' ...
0115             'detected: using average of the z values.']);
0116   fs = find(ind(1:end-1) == 0 & ind(2:end) == 1);
0117   fe = find(ind(1:end-1) == 1 & ind(2:end) == 0);
0118   for i = 1 : length(fs)
0119     % averaging z values
0120     z(fe(i)) = mean(z(fs(i):fe(i)));
0121   end
0122   x = x(~ind(2:end));
0123   y = y(~ind(2:end));
0124   z = z(~ind(2:end));
0125 end
0126
0127 if numel(x) < 3
0128   error('MATLAB:griddata:NotEnoughSamplePts',...
0129         'Not enough unique sample points specified.');
0130 end
0131
0132 switch lower(method),
0133   case 'linear'
0134     zi = linear(x,y,z,xi,yi);
0135   case 'cubic'
0136     zi = cubic(x,y,z,xi,yi);
0137   case 'nearest'
0138     zi = nearest(x,y,z,xi,yi);
0139   case {'invdist','v4'}
0140     zi = gdatav4(x,y,z,xi,yi);
0141   otherwise
0142     error('MATLAB:griddata:UnknownMethod', 'Unknown method.');
0143 end
0144
0145 if nargout<=1, xi = zi; end
0146
0147
0148 %------------------------------------------------------------
0149 function zi = linear(x,y,z,xi,yi)
0150 %LINEAR Triangle-based linear interpolation
0151
0152 %   Reference: David F. Watson, "Contouring: A guide
0153 %   to the analysis and display of spacial data", Pergamon, 1994.
0154
0155
0156 siz = size(xi);
0157 xi = xi(:); yi = yi(:); % Treat these as columns
0158 x = x(:); y = y(:); z = z(:);
0159
0160 dt = DelaunayTri(x,y);
0161 scopedWarnOff = warning('off', 'MATLAB:TriRep:EmptyTri2DWarnId');
0162 restoreWarnOff = onCleanup(@()warning(scopedWarnOff));
0163 dtt = dt.Triangulation;
0164 if isempty(dtt)
0165   error('MATLAB:griddata:EmptyTriangulation','Error computing Delaunay triangulation. The sample datapoints may be collinear.');
0166 end
0167
0168
0169 if(isreal(z))
0170     F = TriScatteredInterp(dt,z);
0171     zi = F(xi,yi);
0172 else
0173     zre = real(z);
0174     zim = imag(z);
0175     F = TriScatteredInterp(dt,zre);
0176     zire = F(xi,yi);
0177     F.V = zim;
0178     ziim = F(xi,yi);
0179     zi = complex(zire,ziim);
0180 end
0181 zi = reshape(zi,siz);
0182
0183
0184
0185
0186 %------------------------------------------------------------
0187
0188 %------------------------------------------------------------
0189 function zi = cubic(x,y,z,xi,yi)
0190 %TRIANGLE Triangle-based cubic interpolation
0191
0192 %   Reference: T. Y. Yang, "Finite Element Structural Analysis",
0193 %   Prentice Hall, 1986.  pp. 446-449.
0194 %
0195 %   Reference: David F. Watson, "Contouring: A guide
0196 %   to the analysis and display of spacial data", Pergamon, 1994.
0197
0198 % Triangularize the data
0199
0200 dt = DelaunayTri([x(:) y(:)]);
0201 scopedWarnOff = warning('off', 'MATLAB:TriRep:EmptyTri2DWarnId');
0202 restoreWarnOff = onCleanup(@()warning(scopedWarnOff));
0203 tri = dt.Triangulation;
0204 if isempty(tri),
0205   error('MATLAB:griddata:EmptyTriangulation','Error computing Delaunay triangulation. The sample datapoints may be collinear.');
0206 end
0207
0208 % Find the enclosing triangle (t)
0209 siz = size(xi);
0210 t = dt.pointLocation(xi(:),yi(:));
0211 t = reshape(t,siz);
0212
0213 if(isreal(z))
0214     zi = cubicmx(x,y,z,xi,yi,tri,t);
0215 else
0216     zre = real(z);
0217     zim = imag(z);
0218     zire = cubicmx(x,y,zre,xi,yi,tri,t);
0219     ziim = cubicmx(x,y,zim,xi,yi,tri,t);
0220     zi = complex(zire,ziim);
0221 end
0222 %------------------------------------------------------------
0223
0224 %------------------------------------------------------------
0225 function zi = nearest(x,y,z,xi,yi)
0226 %NEAREST Triangle-based nearest neightbor interpolation
0227
0228 %   Reference: David F. Watson, "Contouring: A guide
0229 %   to the analysis and display of spacial data", Pergamon, 1994.
0230
0231 siz = size(xi);
0232 xi = xi(:); yi = yi(:); % Treat these a columns
0233 dt = DelaunayTri(x,y);
0234 scopedWarnOff = warning('off', 'MATLAB:TriRep:EmptyTri2DWarnId');
0235 restoreWarnOff = onCleanup(@()warning(scopedWarnOff));
0236 dtt = dt.Triangulation;
0237 if isempty(dtt)
0238   error('MATLAB:griddata:EmptyTriangulation','Error computing Delaunay triangulation. The sample datapoints may be collinear.');
0239 end
0240
0241 k = dt.nearestNeighbor(xi,yi);
0242 zi = k;
0243 d = find(isfinite(k));
0244 zi(d) = z(k(d));
0245 zi = reshape(zi,siz);
0246
0247
0248 %----------------------------------------------------------
0249
0250
0251 %----------------------------------------------------------
0252 function [xi,yi,zi] = gdatav4(x,y,z,xi,yi)
0253 %GDATAV4 MATLAB 4 GRIDDATA interpolation
0254
0255 %   Reference:  David T. Sandwell, Biharmonic spline
0256 %   interpolation of GEOS-3 and SEASAT altimeter
0257 %   data, Geophysical Research Letters, 2, 139-142,
0258 %   1987.  Describes interpolation using value or
0259 %   gradient of value in any dimension.
0260
0261 xy = x(:) + y(:)*sqrt(-1);
0262
0263 % Determine distances between points
0264 d = xy(:,ones(1,length(xy)));
0265 d = abs(d - d.');
0266 n = size(d,1);
0267 % Replace zeros along diagonal with ones (so these don't show up in the
0268 % find below or in the Green's function calculation).
0269 d(1:n+1:numel(d)) = ones(1,n);
0270
0271 non = find(d == 0);
0272 if ~isempty(non),
0273   % If we've made it to here, then some points aren't distinct.  Remove
0274   % the non-distinct points by averaging.
0275   [r,c] = find(d == 0);
0276   k = find(r < c);
0277   r = r(k); c = c(k); % Extract unique (row,col) pairs
0278   v = (z(r) + z(c))/2; % Average non-distinct pairs
0279
0280   rep = find(diff(c)==0);
0281   if ~isempty(rep), % More than two points need to be averaged.
0282     runs = find(diff(diff(c)==0)==1)+1;
0283     for i=1:length(runs),
0284       k = find(c==c(runs(i))); % All the points in a run
0285       v(runs(i)) = mean(z([r(k);c(runs(i))])); % Average (again)
0286     end
0287   end
0288   z(r) = v;
0289   if ~isempty(rep),
0290     z(r(runs)) = v(runs); % Make sure average is in the dataset
0291   end
0292
0293   % Now remove the extra points.
0294   x(c) = [];
0295   y(c) = [];
0296   z(c) = [];
0297   xy(c,:) = [];
0298   xy(:,c) = [];
0299   d(c,:) = [];
0300   d(:,c) = [];
0301
0302   % Determine the non distinct points
0303   ndp = sort([r;c]);
0304   ndp(find(ndp(1:length(ndp)-1)==ndp(2:length(ndp)))) = [];
0305
0306   warning('MATLAB:griddata:NonDistinctPoints',['Averaged %d non-distinct ' ...
0307             'points.\n         Indices are: %s.'],length(ndp),num2str(ndp'))
0308 end
0309
0310 % Determine weights for interpolation
0311 g = (d.^2) .* (log(d)-1);   % Green's function.
0312 % Fixup value of Green's function along diagonal
0313 g(1:size(d,1)+1:numel(d)) = zeros(size(d,1),1);
0314 weights = g \ z(:);
0315
0316 [m,n] = size(xi);
0317 zi = zeros(size(xi));
0318 jay = sqrt(-1);
0319 xy = xy.';
0320
0321 % Evaluate at requested points (xi,yi).  Loop to save memory.
0322 for i=1:m
0323   for j=1:n
0324     d = abs(xi(i,j)+jay*yi(i,j) - xy);
0325     mask = find(d == 0);
0326     if length(mask)>0, d(mask) = ones(length(mask),1); end
0327     g = (d.^2) .* (log(d)-1);   % Green's function.
0328     % Value of Green's function at zero
0329     if length(mask)>0, g(mask) = zeros(length(mask),1); end
0330     zi(i,j) = g * weights;
0331   end
0332 end
0333
0334 if nargout<=1,
0335   xi = zi;
0336 end
0337
0338
0339
0340
0341```

Generated on Tue 10-May-2016 16:37:45 by m2html © 2005