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