matlab实现三维表面拟合griffit参数选择

输入参数

x,y,z

x，y，z分别为等长度向量，且大于3，对x和y的唯一约束是它们不能都落在x-y平面上的一条直线上

向量长度不等时，报错，输出：'Data vectors are incompatible in size.'

向量长度<3 ，报错，输出：'Insufficient data for surface estimation.'

向量中含有nan时，自动舍去。

 xnodes

数据类型：可以为向量，也可为一个标量整数，若为一个整数，表示x最大值与最小值之间等间距节点的数量。

定义网格中x方向的节点。向量不需要等距。xnode中的数不能重复，必须完全跨越数据

xnode中含重复数据，报错，输出：'xnodes and ynodes must be monotone increasing'

如果xnode没有跨域整个x数据，结合参数extend

例：xnodes=1：0.1：1

表示：x方向划分11个网格，网格间距为0.1

 ynodes

与 xnodes类似。

smoothness

决定估计表面的最终平滑度。这里的值越大意味着表面会更光滑。平滑度必须是非负实数。

缺省值：1

参数：如果这个参数是一个长度为2的向量，那么它定义了与x和y变量相关的相对平滑。允许在x维度上应用与y维度不同的平滑量。

interp

插值数据的插值方案

缺省值：triangle——将网格中的每个单元格分割成一个三角形，然后在每个三角形内部进行线性插值

参数：bilinear——双线性插值

           nearest——最近邻插值。一般不选择，但是我将它作为完整性的一个选项包括进来

regularizer

正则化，表数据保真性。我们可以把表面看作一个板，其中板的弯曲刚度由用户指定为一个相对于数据保真度的重要性的数字。一个更硬的板会导致一个更光滑的表面，但适合数据不是很好。我们也可以把正则化器看作是一个扩散问题，其中相对导热系数是提供的。在这里，插值被看作是在给定一组保持在固定温度下的点的情况下，在一个物体中寻找稳定温度分布的问题。外推是线性的。这两种模式都适合于拉普拉斯正则化。

缺省值：gradient——尝试确保所有地方的坡度都尽可能平滑。它与“扩散”选项有细微的不同，在这里，方向导数是有偏差的，以在网格的细胞边界上是平滑的。

参数：diffusion或者laplacian——用有限差分逼近拉普拉斯算子

讯享网

          springs：使用一个spring模型将节点彼此连接，并将数据点连接到网格中的节点。这个选择将使任何外推都尽可能保持不变。倾向于将表面拖向所有数据的平均值，所以系数不宜过大，导致过于平滑

solver

缺省值：normal——使用matlab的反斜杠运算符求解稀疏系统

参数：symmlq——使用matlab的迭代symmlq求解器

           lsqr——使用matlab的迭代LSQR求解器

maxiter

仅适用于迭代求解器-定义迭代求解器的最大迭代次数

缺省值：min(10000,length(xnodes)*length(ynodes)) 

extend

对 xnodes和 ynodes处理，使得其边界为x，y边界

缺省值：warning——给出警告，自动调整 xnodes和 ynodes

参数：never——报错，停止程序

           always——不给出警告，自动调整 xnodes和 ynodes

tilesize

如果网格太大，无法在一个单一的估计步骤中解决，可以将其构建为一组贴图，每次使用更小的网格面积

如果数据非常稀疏，以至于一些贴图中不包含足够的数据来建模，那么这些贴图将被保留为nan。

缺省值：inf

参数：建议为100-200，最小值为3

overlap

缺省值：0.2

参数：         0 <= overlap <= 0.5

autoscale

一些数据在各自的x轴和y轴上可能有很大不同的尺度。如果发生这种情况，那么规范化可能会遇到困难

缺省值：on——将导致gridfit将x和y节点间隔缩放为单位长度

参数：off——禁用自动缩放

代码

function [zgrid,xgrid,ygrid] = gridfit(x,y,z,xnodes,ynodes,varargin) % Copyright 2018 The MathWorks, Inc. % 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. % x,y,z 长度一致， % ANY points containing a NaN are ignored in the estimation % NAN会自动忽略 % % 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 or vector of length 2 - 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. % % If this parameter is a vector of length 2, then it defines % the relative smoothing to be associated with the x and y % variables. This allows the user to apply a different amount % of smoothing in the x dimension compared to the y dimension. % % 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: The 'springs' regularizer tends to drag the surface % towards the mean of all the data, so too large a smoothing % parameter may be a problem. % % % '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:  % Release: 2.0 % Release date: 5/23/06 % set defaults params.smoothness = 1; params.interp = 'triangle'; params.regularizer = 'gradient'; params.solver = 'backslash'; 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); end end % check the parameters for acceptability params = check_params(params); % ensure all of x,y,z,xnodes,ynodes are column vectors, % also drop any NaN data x=x(:); y=y(:); z=z(:); k = isnan(x) | isnan(y) | isnan(z); if any(k) x(k)=[]; y(k)=[]; z(k)=[]; end xmin = min(x); xmax = max(x); ymin = min(y); ymax = max(y); % did they supply a scalar for the nodes? if length(xnodes)==1 xnodes = linspace(xmin,xmax,xnodes)'; xnodes(end) = xmax; % make sure it hits the max end if length(ynodes)==1 ynodes = linspace(ymin,ymax,ynodes)'; ynodes(end) = ymax; % make sure it hits the max end xnodes=xnodes(:); ynodes=ynodes(:); dx = diff(xnodes); dy = diff(ynodes); nx = length(xnodes); ny = length(ynodes); ngrid = nx*ny; % set the scaling if autoscale was on if strcmpi(params.autoscale,'on') params.xscale = mean(dx); params.yscale = mean(dy); params.autoscale = 'off'; end % check to see if any tiling is necessary if (params.tilesize < max(nx,ny)) % split it into smaller tiles. compute zgrid and ygrid % at the very end if requested zgrid = tiled_gridfit(x,y,z,xnodes,ynodes,params); else % its a single tile. % mask must be either an empty array, or a boolean % aray of the same size as the final grid. nmask = size(params.mask); if ~isempty(params.mask) && ((nmask(2)~=nx) || (nmask(1)~=ny)) if ((nmask(2)==ny) || (nmask(1)==nx)) error 'Mask array is probably transposed from proper orientation.' else error 'Mask array must be the same size as the final grid.' end end if ~isempty(params.mask) params.maskflag = 1; else params.maskflag = 0; end % default for maxiter? if isempty(params.maxiter) params.maxiter = min(10000,nx*ny); end % check lengths of the data n = length(x); if (length(y)~=n) || (length(z)~=n) error 'Data vectors are incompatible in size.' end if n<3 error 'Insufficient data for surface estimation.' end % verify the nodes are distinct if any(diff(xnodes)<=0) || any(diff(ynodes)<=0) error 'xnodes and ynodes must be monotone increasing' end % do we need to tweak the first or last node in x or y? if xmin<xnodes(1) switch params.extend case 'always' xnodes(1) = xmin; case 'warning' warning('GRIDFIT:extend',['xnodes(1) was decreased by: ',num2str(xnodes(1)-xmin),', new node = ',num2str(xmin)]) xnodes(1) = xmin; case 'never' error(['Some x (',num2str(xmin),') falls below xnodes(1) by: ',num2str(xnodes(1)-xmin)]) end end if xmax>xnodes(end) switch params.extend case 'always' xnodes(end) = xmax; case 'warning' warning('GRIDFIT:extend',['xnodes(end) was increased by: ',num2str(xmax-xnodes(end)),', new node = ',num2str(xmax)]) xnodes(end) = xmax; case 'never' error(['Some x (',num2str(xmax),') falls above xnodes(end) by: ',num2str(xmax-xnodes(end))]) end end if ymin<ynodes(1) switch params.extend case 'always' ynodes(1) = ymin; case 'warning' warning('GRIDFIT:extend',['ynodes(1) was decreased by: ',num2str(ynodes(1)-ymin),', new node = ',num2str(ymin)]) ynodes(1) = ymin; case 'never' error(['Some y (',num2str(ymin),') falls below ynodes(1) by: ',num2str(ynodes(1)-ymin)]) end end if ymax>ynodes(end) switch params.extend case 'always' ynodes(end) = ymax; case 'warning' warning('GRIDFIT:extend',['ynodes(end) was increased by: ',num2str(ymax-ynodes(end)),', new node = ',num2str(ymax)]) ynodes(end) = ymax; case 'never' error(['Some y (',num2str(ymax),') falls above ynodes(end) by: ',num2str(ymax-ynodes(end))]) end end % determine which cell in the array each point lies in [junk,indx] = histc(x,xnodes); %#ok [junk,indy] = histc(y,ynodes); %#ok % any point falling at the last node is taken to be % inside the last cell in x or y. k=(indx==nx); indx(k)=indx(k)-1; k=(indy==ny); indy(k)=indy(k)-1; ind = indy + ny*(indx-1); % Do we have a mask to apply? if params.maskflag % if we do, then we need to ensure that every % cell with at least one data point also has at % least all of its corners unmasked. params.mask(ind) = 1; params.mask(ind+1) = 1; params.mask(ind+ny) = 1; params.mask(ind+ny+1) = 1; end % interpolation equations for each point tx = min(1,max(0,(x - xnodes(indx))./dx(indx))); ty = min(1,max(0,(y - ynodes(indy))./dy(indy))); % Future enhancement: add cubic interpolant switch params.interp case 'triangle' % linear interpolation inside each triangle k = (tx > ty); L = ones(n,1); L(k) = ny; t1 = min(tx,ty); t2 = max(tx,ty); A = sparse(repmat((1:n)',1,3),[ind,ind+ny+1,ind+L], ... [1-t2,t1,t2-t1],n,ngrid); case 'nearest' % nearest neighbor interpolation in a cell k = round(1-ty) + round(1-tx)*ny; A = sparse((1:n)',ind+k,ones(n,1),n,ngrid); case 'bilinear' % bilinear interpolation in a cell A = sparse(repmat((1:n)',1,4),[ind,ind+1,ind+ny,ind+ny+1], ... [(1-tx).*(1-ty), (1-tx).*ty, tx.*(1-ty), tx.*ty], ... n,ngrid); end rhs = z; % do we have relative smoothing parameters? if numel(params.smoothness) == 1 % it was scalar, so treat both dimensions equally smoothparam = params.smoothness; xyRelativeStiffness = [1;1]; else % It was a vector, so anisotropy reigns. % I've already checked that the vector was of length 2 smoothparam = sqrt(prod(params.smoothness)); xyRelativeStiffness = params.smoothness(:)./smoothparam; end % Build regularizer. Add del^4 regularizer one day. switch params.regularizer case 'springs' % zero "rest length" springs [i,j] = meshgrid(1:nx,1:(ny-1)); ind = j(:) + ny*(i(:)-1); m = nx*(ny-1); stiffness = 1./(dy/params.yscale); Areg = sparse(repmat((1:m)',1,2),[ind,ind+1], ... xyRelativeStiffness(2)*stiffness(j(:))*[-1 1], ... m,ngrid); [i,j] = meshgrid(1:(nx-1),1:ny); ind = j(:) + ny*(i(:)-1); m = (nx-1)*ny; stiffness = 1./(dx/params.xscale); Areg = [Areg;sparse(repmat((1:m)',1,2),[ind,ind+ny], ... xyRelativeStiffness(1)*stiffness(i(:))*[-1 1],m,ngrid)]; [i,j] = meshgrid(1:(nx-1),1:(ny-1)); ind = j(:) + ny*(i(:)-1); m = (nx-1)*(ny-1); stiffness = 1./sqrt((dx(i(:))/params.xscale/xyRelativeStiffness(1)).^2 + ... (dy(j(:))/params.yscale/xyRelativeStiffness(2)).^2); Areg = [Areg;sparse(repmat((1:m)',1,2),[ind,ind+ny+1], ... stiffness*[-1 1],m,ngrid)]; Areg = [Areg;sparse(repmat((1:m)',1,2),[ind+1,ind+ny], ... stiffness*[-1 1],m,ngrid)]; case {'diffusion' 'laplacian'} % thermal diffusion using Laplacian (del^2) [i,j] = meshgrid(1:nx,2:(ny-1)); ind = j(:) + ny*(i(:)-1); dy1 = dy(j(:)-1)/params.yscale; dy2 = dy(j(:))/params.yscale; Areg = sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ... xyRelativeStiffness(2)*[-2./(dy1.*(dy1+dy2)), ... 2./(dy1.*dy2), -2./(dy2.*(dy1+dy2))],ngrid,ngrid); [i,j] = meshgrid(2:(nx-1),1:ny); ind = j(:) + ny*(i(:)-1); dx1 = dx(i(:)-1)/params.xscale; dx2 = dx(i(:))/params.xscale; Areg = Areg + sparse(repmat(ind,1,3),[ind-ny,ind,ind+ny], ... xyRelativeStiffness(1)*[-2./(dx1.*(dx1+dx2)), ... 2./(dx1.*dx2), -2./(dx2.*(dx1+dx2))],ngrid,ngrid); case 'gradient' % Subtly different from the Laplacian. A point for future % enhancement is to do it better for the triangle interpolation % case. [i,j] = meshgrid(1:nx,2:(ny-1)); ind = j(:) + ny*(i(:)-1); dy1 = dy(j(:)-1)/params.yscale; dy2 = dy(j(:))/params.yscale; Areg = sparse(repmat(ind,1,3),[ind-1,ind,ind+1], ... xyRelativeStiffness(2)*[-2./(dy1.*(dy1+dy2)), ... 2./(dy1.*dy2), -2./(dy2.*(dy1+dy2))],ngrid,ngrid); [i,j] = meshgrid(2:(nx-1),1:ny); ind = j(:) + ny*(i(:)-1); dx1 = dx(i(:)-1)/params.xscale; dx2 = dx(i(:))/params.xscale; Areg = [Areg;sparse(repmat(ind,1,3),[ind-ny,ind,ind+ny], ... xyRelativeStiffness(1)*[-2./(dx1.*(dx1+dx2)), ... 2./(dx1.*dx2), -2./(dx2.*(dx1+dx2))],ngrid,ngrid)]; end nreg = size(Areg,1); % 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)]; rhs = [rhs;zeros(nreg,1)]; % do we have a mask to apply? if params.maskflag unmasked = find(params.mask); end % solve the full system, with regularizer attached switch params.solver case {'\' 'backslash'} if params.maskflag % there is a mask to use zgrid=nan(ny,nx); zgrid(unmasked) = A(:,unmasked)\rhs; else % no mask zgrid = reshape(A\rhs,ny,nx); end case 'normal' % The normal equations, solved with \. Can be faster % for huge numbers of data points, but reasonably % sized grids. The regularizer makes A well conditioned % so the normal equations are not a terribly bad thing % here. if params.maskflag % there is a mask to use Aunmasked = A(:,unmasked); zgrid=nan(ny,nx); zgrid(unmasked) = (Aunmasked'*Aunmasked)\(Aunmasked'*rhs); else zgrid = reshape((A'*A)\(A'*rhs),ny,nx); end case 'symmlq' % iterative solver - symmlq - requires a symmetric matrix, % so use it to solve the normal equations. No preconditioner. tol = abs(max(z)-min(z))*1.e-13; if params.maskflag % there is a mask to use zgrid=nan(ny,nx); [zgrid(unmasked),flag] = symmlq(A(:,unmasked)'*A(:,unmasked), ... A(:,unmasked)'*rhs,tol,params.maxiter); else [zgrid,flag] = symmlq(A'*A,A'*rhs,tol,params.maxiter); zgrid = reshape(zgrid,ny,nx); end % display a warning if convergence problems switch flag case 0 % no problems with convergence case 1 % SYMMLQ iterated MAXIT times but did not converge. warning('GRIDFIT:solver',['Symmlq performed ',num2str(params.maxiter), ... ' iterations but did not converge.']) case 3 % SYMMLQ stagnated, successive iterates were the same warning('GRIDFIT:solver','Symmlq stagnated without apparent convergence.') otherwise warning('GRIDFIT:solver',['One of the scalar quantities calculated in',... ' symmlq was too small or too large to continue computing.']) end case 'lsqr' % iterative solver - lsqr. No preconditioner here. tol = abs(max(z)-min(z))*1.e-13; if params.maskflag % there is a mask to use zgrid=nan(ny,nx); [zgrid(unmasked),flag] = lsqr(A(:,unmasked),rhs,tol,params.maxiter); else [zgrid,flag] = lsqr(A,rhs,tol,params.maxiter); zgrid = reshape(zgrid,ny,nx); end % display a warning if convergence problems switch flag case 0 % no problems with convergence case 1 % lsqr iterated MAXIT times but did not converge. warning('GRIDFIT:solver',['Lsqr performed ', ... num2str(params.maxiter),' iterations but did not converge.']) case 3 % lsqr stagnated, successive iterates were the same warning('GRIDFIT:solver','Lsqr stagnated without apparent convergence.') case 4 warning('GRIDFIT:solver',['One of the scalar quantities calculated in',... ' LSQR was too small or too large to continue computing.']) end end % switch params.solver end % if params.tilesize... % only generate xgrid and ygrid if requested. if nargout>1 [xgrid,ygrid]=meshgrid(xnodes,ynodes); end % ============================================ % End of main function - gridfit % ============================================ % ============================================ % subfunction - parse_pv_pairs % ============================================ function params=parse_pv_pairs(params,pv_pairs) % parse_pv_pairs: parses sets of property value pairs, allows defaults % usage: params=parse_pv_pairs(default_params,pv_pairs) % % arguments: (input) % default_params - structure, with one field for every potential % property/value pair. Each field will contain the default % value for that property. If no default is supplied for a % given property, then that field must be empty. % % pv_array - cell array of property/value pairs. % Case is ignored when comparing properties to the list % of field names. Also, any unambiguous shortening of a % field/property name is allowed. % % arguments: (output) % params - parameter struct that reflects any updated property/value % pairs in the pv_array. % % Example usage: % First, set default values for the parameters. Assume we % have four parameters that we wish to use optionally in % the function examplefun. % % - 'viscosity', which will have a default value of 1 % - 'volume', which will default to 1 % - 'pie' - which will have default value 3.9793 % - 'description' - a text field, left empty by default % % The first argument to examplefun is one which will always be % supplied. % % function examplefun(dummyarg1,varargin) % params.Viscosity = 1; % params.Volume = 1; % params.Pie = 3.9793 % % params.Description = ''; % params=parse_pv_pairs(params,varargin); % params % % Use examplefun, overriding the defaults for 'pie', 'viscosity' % and 'description'. The 'volume' parameter is left at its default. % % examplefun(rand(10),'vis',10,'pie',3,'Description','Hello world') % % params = % Viscosity: 10 % Volume: 1 % Pie: 3 % Description: 'Hello world' % % Note that capitalization was ignored, and the property 'viscosity' % was truncated as supplied. Also note that the order the pairs were % supplied was arbitrary. npv = length(pv_pairs); n = npv/2; if n~=floor(n) error 'Property/value pairs must come in PAIRS.' end if n<=0 % just return the defaults return end if ~isstruct(params) error 'No structure for defaults was supplied' end % there was at least one pv pair. process any supplied propnames = fieldnames(params); lpropnames = lower(propnames); for i=1:n p_i = lower(pv_pairs{2*i-1}); v_i = pv_pairs{2*i}; ind = strmatch(p_i,lpropnames,'exact'); if isempty(ind) ind = find(strncmp(p_i,lpropnames,length(p_i))); if isempty(ind) error(['No matching property found for: ',pv_pairs{2*i-1}]) elseif length(ind)>1 error(['Ambiguous property name: ',pv_pairs{2*i-1}]) end end p_i = propnames{ind}; % override the corresponding default in params params = setfield(params,p_i,v_i); %#ok end % ============================================ % subfunction - check_params % ============================================ function params = check_params(params) % check the parameters for acceptability % smoothness == 1 by default if isempty(params.smoothness) params.smoothness = 1; else if (numel(params.smoothness)>2) || any(params.smoothness<=0) error 'Smoothness must be scalar (or length 2 vector), real, finite, and positive.' end end % regularizer - must be one of 4 options - the second and % third are actually synonyms. valid = {'springs', 'diffusion', 'laplacian', 'gradient'}; if isempty(params.regularizer) params.regularizer = 'diffusion'; end ind = find(strncmpi(params.regularizer,valid,length(params.regularizer))); if (length(ind)==1) params.regularizer = valid{ind}; else error(['Invalid regularization method: ',params.regularizer]) end % interp must be one of: % 'bilinear', 'nearest', or 'triangle' % but accept any shortening thereof. valid = {'bilinear', 'nearest', 'triangle'}; if isempty(params.interp) params.interp = 'triangle'; end ind = find(strncmpi(params.interp,valid,length(params.interp))); if (length(ind)==1) params.interp = valid{ind}; else error(['Invalid interpolation method: ',params.interp]) end % solver must be one of: % 'backslash', '\', 'symmlq', 'lsqr', or 'normal' % but accept any shortening thereof. valid = {'backslash', '\', 'symmlq', 'lsqr', 'normal'}; if isempty(params.solver) params.solver = '\'; end ind = find(strncmpi(params.solver,valid,length(params.solver))); if (length(ind)==1) params.solver = valid{ind}; else error(['Invalid solver option: ',params.solver]) end % extend must be one of: % 'never', 'warning', 'always' % but accept any shortening thereof. valid = {'never', 'warning', 'always'}; if isempty(params.extend) params.extend = 'warning'; end ind = find(strncmpi(params.extend,valid,length(params.extend))); if (length(ind)==1) params.extend = valid{ind}; else error(['Invalid extend option: ',params.extend]) end % tilesize == inf by default if isempty(params.tilesize) params.tilesize = inf; elseif (length(params.tilesize)>1) || (params.tilesize<3) error 'Tilesize must be scalar and > 0.' end % overlap == 0.20 by default if isempty(params.overlap) params.overlap = 0.20; elseif (length(params.overlap)>1) || (params.overlap<0) || (params.overlap>0.5) error 'Overlap must be scalar and 0 < overlap < 1.' end % ============================================ % subfunction - tiled_gridfit % ============================================ function zgrid=tiled_gridfit(x,y,z,xnodes,ynodes,params) % tiled_gridfit: a tiled version of gridfit, continuous across tile boundaries % usage: [zgrid,xgrid,ygrid]=tiled_gridfit(x,y,z,xnodes,ynodes,params) % % Tiled_gridfit is used when the total grid is far too large % to model using a single call to gridfit. While gridfit may take % only a second or so to build a 100x100 grid, a 2000x2000 grid % will probably not run at all due to memory problems. % % Tiles in the grid with insufficient data (<4 points) will be % filled with NaNs. Avoid use of too small tiles, especially % if your data has holes in it that may encompass an entire tile. % % A mask may also be applied, in which case tiled_gridfit will % subdivide the mask into tiles. Note that any boolean mask % provided is assumed to be the size of the complete grid. % % Tiled_gridfit may not be fast on huge grids, but it should run % as long as you use a reasonable tilesize. 8-) % Note that we have already verified all parameters in check_params % Matrix elements in a square tile tilesize = params.tilesize; % Size of overlap in terms of matrix elements. Overlaps % of purely zero cause problems, so force at least two % elements to overlap. overlap = max(2,floor(tilesize*params.overlap)); % reset the tilesize for each particular tile to be inf, so % we will never see a recursive call to tiled_gridfit Tparams = params; Tparams.tilesize = inf; nx = length(xnodes); ny = length(ynodes); zgrid = zeros(ny,nx); % linear ramp for the bilinear interpolation rampfun = inline('(t-t(1))/(t(end)-t(1))','t'); % loop over each tile in the grid h = waitbar(0,'Relax and have a cup of JAVA. Its my treat.'); warncount = 0; xtind = 1:min(nx,tilesize); while ~isempty(xtind) && (xtind(1)<=nx) xinterp = ones(1,length(xtind)); if (xtind(1) ~= 1) xinterp(1:overlap) = rampfun(xnodes(xtind(1:overlap))); end if (xtind(end) ~= nx) xinterp((end-overlap+1):end) = 1-rampfun(xnodes(xtind((end-overlap+1):end))); end ytind = 1:min(ny,tilesize); while ~isempty(ytind) && (ytind(1)<=ny) % update the waitbar waitbar((xtind(end)-tilesize)/nx + tilesize*ytind(end)/ny/nx) yinterp = ones(length(ytind),1); if (ytind(1) ~= 1) yinterp(1:overlap) = rampfun(ynodes(ytind(1:overlap))); end if (ytind(end) ~= ny) yinterp((end-overlap+1):end) = 1-rampfun(ynodes(ytind((end-overlap+1):end))); end % was a mask supplied? if ~isempty(params.mask) submask = params.mask(ytind,xtind); Tparams.mask = submask; end % extract data that lies in this grid tile k = (x>=xnodes(xtind(1))) & (x<=xnodes(xtind(end))) & ... (y>=ynodes(ytind(1))) & (y<=ynodes(ytind(end))); k = find(k); if length(k)<4 if warncount == 0 warning('GRIDFIT:tiling','A tile was too underpopulated to model. Filled with NaNs.') end warncount = warncount + 1; % fill this part of the grid with NaNs zgrid(ytind,xtind) = NaN; else % build this tile zgtile = gridfit(x(k),y(k),z(k),xnodes(xtind),ynodes(ytind),Tparams); % bilinear interpolation (using an outer product) interp_coef = yinterp*xinterp; % accumulate the tile into the complete grid zgrid(ytind,xtind) = zgrid(ytind,xtind) + zgtile.*interp_coef; end % step to the next tile in y if ytind(end)<ny ytind = ytind + tilesize - overlap; % are we within overlap elements of the edge of the grid? if (ytind(end)+max(3,overlap))>=ny % extend this tile to the edge ytind = ytind(1):ny; end else ytind = ny+1; end end % while loop over y % step to the next tile in x if xtind(end)<nx xtind = xtind + tilesize - overlap; % are we within overlap elements of the edge of the grid? if (xtind(end)+max(3,overlap))>=nx % extend this tile to the edge xtind = xtind(1):nx; end else xtind = nx+1; end end % while loop over x % close down the waitbar close(h) if warncount>0 warning('GRIDFIT:tiling',[num2str(warncount),' tiles were underpopulated & filled with NaNs']) end