Mercurial > octave-nkf
view scripts/geometry/griddata.m @ 6823:9fddcc586065
[project @ 2007-08-24 08:27:27 by dbateman]
author | dbateman |
---|---|
date | Fri, 24 Aug 2007 08:27:29 +0000 |
parents | |
children | 8618f29520c6 |
line wrap: on
line source
## Copyright (C) 1999,2000 Kai Habel ## ## This file is part of Octave. ## ## Octave is free software; you can redistribute it and/or modify it ## under the terms of the GNU General Public License as published by ## the Free Software Foundation; either version 2, or (at your option) ## any later version. ## ## Octave is distributed in the hope that it will be useful, but ## WITHOUT ANY WARRANTY; without even the implied warranty of ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ## General Public License for more details. ## ## You should have received a copy of the GNU General Public License ## along with Octave; see the file COPYING. If not, write to the Free ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA ## 02110-1301, USA. ## -*- texinfo -*- ## @deftypefn {Function File} {@var{zi} =} griddata (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{method}) ## @deftypefnx {Function File} {[@var{xi}, @var{yi}, @var{zi}] =} griddata (@var{x}, @var{y}, @var{z}, @var{xi}, @var{yi}, @var{method}) ## ## Generate a regular mesh from irregular data using interpolation. ## The function is defined by @code{@var{z} = f (@var{x}, @var{y})}. ## The interpolation points are all @code{(@var{xi}, @var{yi})}. If ## @var{xi}, @var{yi} are vectors then they are made into a 2D mesh. ## ## The interpolation method can be 'nearest', 'cubic' or 'linear'. ## If method is omitted it defaults to 'linear'. ## @seealso{delaunay} ## @end deftypefn ## Author: Kai Habel <kai.habel@gmx.de> ## Adapted-by: Alexander Barth <barth.alexander@gmail.com> ## xi and yi are not "meshgridded" if both are vectors ## of the same size (for compatibility) function [rx, ry, rz] = griddata (x,y,z,xi,yi,method) if (nargin == 5) method="linear"; endif if (nargin < 5 || nargin > 7) print_usage(); endif if (ischar(method)) method=tolower(method); endif if (!all( (size(x)==size(y)) & (size(x)==size(z)))) error("griddata: x,y,z must be vectors of same length"); endif ## meshgrid xi and yi if they are vectors unless they ## are vectors of the same length if (isvector(xi) && isvector(yi) && numel(xi) ~= numel(yi)) [xi,yi]=meshgrid(xi,yi); endif if (any(size(xi) != size(yi))) error("griddata: xi and yi must be vectors or matrices of same size"); endif [nr,nc] = size(xi); ## triangulate data tri = delaunay(x,y); zi = nan(size(xi)); if strcmp(method,"cubic") error("griddata(...,'cubic') cubic interpolation not yet implemented\n") elseif strcmp(method,'nearest') ## search index of nearest point idx = dsearch(x,y,tri,xi,yi); valid = !isnan(idx); zi(valid) = z(idx(valid)); elseif strcmp(method,'linear') ## search for every point the enclosing triangle tri_list= tsearch(x,y,tri,xi(:),yi(:)); ## only keep the points within triangles. valid = !isnan(reshape(tri_list,size(xi))); tri_list = tri_list(!isnan(tri_list)); nr_t = rows(tri_list); ## assign x,y,z for each point of triangle x1 = x(tri(tri_list,1));y1=y(tri(tri_list,1));z1=z(tri(tri_list,1)); x2 = x(tri(tri_list,2));y2=y(tri(tri_list,2));z2=z(tri(tri_list,2)); x3 = x(tri(tri_list,3));y3=y(tri(tri_list,3));z3=z(tri(tri_list,3)); ## calculate norm vector N = cross([x2-x1, y2-y1, z2-z1],[x3-x1, y3-y1, z3-z1]); N_norm = sqrt(sumsq(N,2)); N = N ./ N_norm(:,[1,1,1]); ## calculate D of plane equation ## Ax+By+Cz+D=0; D = -(N(:,1) .* x1 + N(:,2) .* y1 + N(:,3) .* z1); ## calculate zi by solving plane equation for xi,yi zi(valid) = -(N(:,1).*xi(valid) + N(:,2).*yi(valid) + D ) ./ N(:,3); else error("griddata: unknown interpolation method"); endif if nargout == 3 rx = xi; ry = yi; rz = zi; elseif nargout == 1 rx = zi; elseif nargout == 0 mesh(xi, yi, zi); endif endfunction %!test %! [xx,yy]=meshgrid(linspace(-1,1,32)); %! x = xx(:); %! x = x + 10 * (2 * round(rand(size(x))) - 1) * eps; %! y = yy(:); %! y = y + 10 * (2 * round(rand(size(y))) - 1) * eps; %! z = sin(2*(x.^2+y.^2)); %! zz = griddata(x,y,z,xx,yy,'linear'); %! zz2 = sin(2*(xx.^2+yy.^2)); %! zz2(isnan(zz)) = NaN; %! assert (zz, zz2, 100 * eps) %!demo %! x=2*rand(100,1)-1; %! y=2*rand(size(x))-1; %! z=sin(2*(x.^2+y.^2)); %! [xx,yy]=meshgrid(linspace(-1,1,32)); %! griddata(x,y,z,xx,yy); %! title('nonuniform grid sampled at 100 points'); %!demo %! x=2*rand(1000,1)-1; %! y=2*rand(size(x))-1; %! z=sin(2*(x.^2+y.^2)); %! [xx,yy]=meshgrid(linspace(-1,1,32)); %! griddata(x,y,z,xx,yy); %! title('nonuniform grid sampled at 1000 points'); %!demo %! x=2*rand(1000,1)-1; %! y=2*rand(size(x))-1; %! z=sin(2*(x.^2+y.^2)); %! [xx,yy]=meshgrid(linspace(-1,1,32)); %! griddata(x,y,z,xx,yy,'nearest'); %! title('nonuniform grid sampled at 1000 points with nearest neighbor');