Mercurial > hg > octave-nkf

comparison scripts/general/interpn.m @ 6702:b2391d403ed2

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
[project @ 2007-06-12 21:39:26 by dbateman]
author dbateman
date Tue, 12 Jun 2007 21:39:27 +0000
parents
children 01036667884a
comparison
equal deleted inserted replaced
6701:3933e0693fe0 6702:b2391d403ed2
1 ## Copyright (C) 2007 David Bateman
2 ##
3 ## This file is part of Octave.
4 ##
5 ## Octave is free software; you can redistribute it and/or modify it
6 ## under the terms of the GNU General Public License as published by
7 ## the Free Software Foundation; either version 2, or (at your option)
8 ## any later version.
9 ##
10 ## Octave is distributed in the hope that it will be useful, but
11 ## WITHOUT ANY WARRANTY; without even the implied warranty of
12 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 ## General Public License for more details.
14 ##
15 ## You should have received a copy of the GNU General Public License
16 ## along with Octave; see the file COPYING. If not, write to the Free
17 ## Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
18 ## 02110-1301, USA.
19
20 ## -*- texinfo -*-
21 ## @deftypefn {Function File} {@var{vi} =} interpn (@var{x1}, @var{x2}, @dots{}, @var{v}, @var{y1}, @var{y2}, @dots{})
22 ## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}, @var{y1}, @var{y2}, @dots{})
23 ## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}, @var{m})
24 ## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v})
25 ## @deftypefnx {Function File} {@var{vi} =} interpn (@dots{}, @var{method})
26 ## @deftypefnx {Function File} {@var{vi} =} interpn (@dots{}, @var{method}, @var{extrapval})
27 ##
28 ## Perform @var{n}-dimensional interpolation, where @var{n} is at least two.
29 ## Each element of then @var{n}-dimensional array @var{v} represents a value
30 ## at a location given by the parameters @var{x1}, @var{x2}, @dots{}, @var{xn}.
31 ## The parameters @var{x1}, @var{x2}, @dots{}, @var{xn} are either
32 ## @var{n}-dimensional arrays of the same size as the array @var{v} in
33 ## the 'ndgrid' format or vectors. The parameters @var{y1}, etc respect a
34 ## similar format to @var{x1}, etc, and they represent the points at which
35 ## the array @var{vi} is interpolated.
36 ##
37 ## If @var{x1}, @dots{}, @var{xn} are ommitted, they are assumed to be
38 ## @code{x1 = 1 : size (@var{v}, 1)}, etc. If @var{m} is specified, then
39 ## the interpolation adds a point half way between each of the interplation
40 ## points. This process is performed @var{m} times. If only @var{v} is
41 ## specified, then @var{m} is assumed to be @code{1}.
42 ##
43 ## Method is one of:
44 ##
45 ## @table @asis
46 ## @item 'nearest'
47 ## Return the nearest neighbour.
48 ## @item 'linear'
49 ## Linear interpolation from nearest neighbours.
50 ## @item 'cubic'
51 ## Cubic interpolation from four nearest neighbours (not implemented yet).
52 ## @item 'spline'
53 ## Cubic spline interpolation--smooth first and second derivatives
54 ## throughout the curve.
55 ## @end table
56 ##
57 ## The default method is 'linear'.
58 ##
59 ## If @var{extrap} is the string 'extrap', then extrapolate values beyond
60 ## the endpoints. If @var{extrap} is a number, replace values beyond the
61 ## endpoints with that number. If @var{extrap} is missing, assume NaN.
62 ## @seealso{interp1, interp2, spline, ndgrid}
63 ## @end deftypefn
64
65 function vi = interpn (varargin)
66
67 method = "linear";
68 extrapval = NaN;
69 nargs = nargin;
70
71 if (nargin < 1)
72 print_usage ();
73 endif
74
75 if (ischar (varargin {end}))
76 method = varargin {end};
77 nargs = nargs - 1;
78 elseif (ischar (varargin {end - 1}))
79 if (! isnumeric (vargin {end}) || ! isscalar (vargin {end}))
80 error ("extrapal is expected to be a numeric scalar");
81 endif
82 method = varargin {end - 1};
83 nargs = nargs - 2;
84 endif
85
86 if (nargs < 3)
87 v = varargin {1};
88 m = 1;
89 if (nargs == 2)
90 m = varargin {2};
91 if (! isnumeric (m) || ! isscalar (m) || floor (m) != m)
92 error ("m is expected to be a integer scalar");
93 endif
94 endif
95 sz = size (v);
96 nd = ndims (v);
97 x = cell (1, nd);
98 y = cell (1, nd);
99 for i = 1 : nd;
100 x{i} = 1 : sz(i);
101 y{i} = 1 : (1 / (2 ^ m)) : sz(i);
102 endfor
103 elseif (! isvector (varargin {1}) && nargs == (ndims (varargin {1}) + 1))
104 v = varargin {1};
105 sz = size (v);
106 nd = ndims (v);
107 x = cell (1, nd);
108 y = varargin (2 : nargs);
109 for i = 1 : nd;
110 x{i} = 1 : sz(i);
111 endfor
112 elseif (rem (nargs, 2) == 1 && nargs ==
113 (2 * ndims (varargin {ceil (nargs / 2)})) + 1)
114 nv = ceil (nargs / 2);
115 v = varargin {nv};
116 sz = size (v);
117 nd = ndims (v);
118 x = varargin (1 : (nv - 1));
119 y = varargin ((nv + 1) : nargs);
120 else
121 error ("wrong number or incorrectly formatted input arguments");
122 endif
123
124 if (any (! cellfun (@isvector, x)))
125 for i = 2 : nd
126 if (! size_equal (x{1}, x{i}) || ! size_equal (x{i}, v))
127 error ("dimensional mismatch");
128 endif
129 idx (1 : nd) = {1};
130 idx (i) = ":";
131 x{i} = x{i}(idx{:});
132 endfor
133 idx (1 : nd) = {1};
134 idx (1) = ":";
135 x{1} = x{1}(idx{:});
136 endif
137
138 if (strcmp (method, "linear") || strcmp (method, "nearest"))
139 if (all (cellfun (@isvector, y)))
140 [y{:}] = ndgrid (y{:});
141 endif
142 elseif (any (! cellfun (@isvector, x)))
143 for i = 1 : nd
144 idx (1 : nd) = {1};
145 idx (i) = ":";
146 y{i} = y{i}(idx{:});
147 endfor
148 endif
149
150 method = tolower (method);
151 if (strcmp (method, "linear"))
152 vi = __lin_interpn__ (x{:}, v, y{:});
153 vi (vi == NaN) = extrapval;
154 elseif (strcmp (method, "nearest"))
155 yshape = size (y{1});
156 yidx = cell (1, nd);
157 for i = 1 : nd
158 y{i} = y{i}(:);
159 yidx{i} = lookup (x{i}(2:end-1), y{i}) + 1;
160 endfor
161 idx = cell (1,nd);
162 for i = 1 : nd
163 idx {i} = yidx{i} + (y{i} - x{i}(yidx{i}).' > ...
164 x{i}(yidx{i} + 1).' - y{i});
165 endfor
166 vi = v (sub2ind (sz, idx{:}));
167 idx = zeros (prod(yshape),1);
168 for i = 1 : nd
169 idx |= y{i} < min (x{i}(:)) | y{i} > max (x{i}(:));
170 endfor
171 vi(idx) = extrapval;
172 vi = reshape (vi, yshape);
173 elseif (strcmp (method, "spline"))
174 vi = __splinen__ (x, v, y, extrapval, "interpn");
175 elseif (strcmp (method, "cubic"))
176 error ("cubic interpolation not yet implemented");
177 else
178 error ("unrecognized interpolation method");
179 endif
180
181 endfunction
182
183 %!demo
184 %! A=[13,-1,12;5,4,3;1,6,2];
185 %! x=[0,1,4]; y=[10,11,12];
186 %! xi=linspace(min(x),max(x),17);
187 %! yi=linspace(min(y),max(y),26)';
188 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"linear").');
189 %! [x,y] = meshgrid(x,y);
190 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
191
192 %!demo
193 %! A=[13,-1,12;5,4,3;1,6,2];
194 %! x=[0,1,4]; y=[10,11,12];
195 %! xi=linspace(min(x),max(x),17);
196 %! yi=linspace(min(y),max(y),26)';
197 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"nearest").');
198 %! [x,y] = meshgrid(x,y);
199 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
200
201 %!#demo
202 %! A=[13,-1,12;5,4,3;1,6,2];
203 %! x=[0,1,2]; y=[10,11,12];
204 %! xi=linspace(min(x),max(x),17);
205 %! yi=linspace(min(y),max(y),26)';
206 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"cubic").');
207 %! [x,y] = meshgrid(x,y);
208 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
209
210 %!demo
211 %! A=[13,-1,12;5,4,3;1,6,2];
212 %! x=[0,1,2]; y=[10,11,12];
213 %! xi=linspace(min(x),max(x),17);
214 %! yi=linspace(min(y),max(y),26)';
215 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"spline").');
216 %! [x,y] = meshgrid(x,y);
217 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off;
218