Mercurial > octave
view scripts/geometry/griddatan.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) 2007 David Bateman <dbateman@free.fr> ## ## 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{yi} =} griddatan (@var{x}, @var{y}, @var{xi}, @var{method}, @var{options}) ## ## Generate a regular mesh from irregular data using interpolation. ## The function is defined by @code{@var{y} = f (@var{x})}. ## The interpolation points are all @var{xi}. ## ## The interpolation method can be 'nearest' or 'linear'. If method is ## omitted it defaults to 'linear'. ## @seealso{griddata, delaunayn} ## @end deftypefn function [yi] = griddatan (x,y,xi,method,varargin) if (nargin == 3) method="linear"; endif if (nargin < 3) print_usage(); endif if (ischar(method)) method=tolower(method); endif [m, n] = size(x); [mi, ni] = size(xi); if (n != ni || size(y,1) != m || size(y,2) != 1) error ("griddatan: dimensional mismatch"); endif ## triangulate data ## tri = delaunayn(x, varargin{:}); tri = delaunayn(x); yi = nan(mi, 1); if strcmp(method,'nearest') ## search index of nearest point idx = dsearchn(x,tri,xi); valid = !isnan(idx); yi(valid) = y(idx(valid)); elseif strcmp(method,'linear') ## search for every point the enclosing triangle [tri_list, bary_list] = tsearchn(x,tri,xi); ## only keep the points within triangles. valid = !isnan(tri_list); tri_list = tri_list(!isnan(tri_list)); bary_list = bary_list(!isnan(tri_list), :); nr_t = rows(tri_list); ## assign x,y for each point of simplex xt = reshape (x(tri(tri_list,:),:),[nr_t, n+1, n]); yt = y(tri(tri_list,:)); ## Use barycentric coordinate of point to calculate yi yi(valid) = sum (y(tri(tri_list,:)) .* bary_list, 2); else method error("griddatan: unknown interpolation method"); endif endfunction %!test %! [xx,yy] = meshgrid(linspace(-1,1,32)); %! xi = [xx(:), yy(:)]; %! x = (2 * rand(100,2) - 1); %! x = [x;1,1;1,-1;-1,-1;-1,1]; %! y = sin(2*(sum(x.^2,2))); %! zz = griddatan(x,y,xi,'linear'); %! zz2 = griddata(x(:,1),x(:,2),y,xi(:,1),xi(:,2),'linear'); %! assert (zz, zz2, 1e-10) %!test %! [xx,yy] = meshgrid(linspace(-1,1,32)); %! xi = [xx(:), yy(:)]; %! x = (2 * rand(100,2) - 1); %! x = [x;1,1;1,-1;-1,-1;-1,1]; %! y = sin(2*(sum(x.^2,2))); %! zz = griddatan(x,y,xi,'nearest'); %! zz2 = griddata(x(:,1),x(:,2),y,xi(:,1),xi(:,2),'nearest'); %! assert (zz, zz2, 1e-10)