Mercurial > hg > octave-lyh
annotate scripts/general/interpn.m @ 8920:eb63fbe60fab
update copyright notices
| author | John W. Eaton <jwe@octave.org> |
|---|---|
| date | Sat, 07 Mar 2009 10:41:27 -0500 |
| parents | fa78cb8d8a5c |
| children | 1bf0ce0930be |
| rev | line source |
|---|---|
| 8920 | 1 ## Copyright (C) 2007, 2008 David Bateman |
| 6702 | 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 | |
| 7016 | 7 ## the Free Software Foundation; either version 3 of the License, or (at |
| 8 ## your option) any later version. | |
| 6702 | 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 | |
| 7016 | 16 ## along with Octave; see the file COPYING. If not, see |
| 17 ## <http://www.gnu.org/licenses/>. | |
| 6702 | 18 |
| 19 ## -*- texinfo -*- | |
| 20 ## @deftypefn {Function File} {@var{vi} =} interpn (@var{x1}, @var{x2}, @dots{}, @var{v}, @var{y1}, @var{y2}, @dots{}) | |
| 21 ## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}, @var{y1}, @var{y2}, @dots{}) | |
| 22 ## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}, @var{m}) | |
| 23 ## @deftypefnx {Function File} {@var{vi} =} interpn (@var{v}) | |
| 24 ## @deftypefnx {Function File} {@var{vi} =} interpn (@dots{}, @var{method}) | |
| 25 ## @deftypefnx {Function File} {@var{vi} =} interpn (@dots{}, @var{method}, @var{extrapval}) | |
| 26 ## | |
| 27 ## Perform @var{n}-dimensional interpolation, where @var{n} is at least two. | |
|
8347
fa78cb8d8a5c
corrections for typos
Brian Gough<bjg@network-theory.co.uk>
parents:
8116
diff
changeset
|
28 ## Each element of the @var{n}-dimensional array @var{v} represents a value |
| 6702 | 29 ## at a location given by the parameters @var{x1}, @var{x2}, @dots{}, @var{xn}. |
| 30 ## The parameters @var{x1}, @var{x2}, @dots{}, @var{xn} are either | |
| 31 ## @var{n}-dimensional arrays of the same size as the array @var{v} in | |
| 32 ## the 'ndgrid' format or vectors. The parameters @var{y1}, etc respect a | |
| 33 ## similar format to @var{x1}, etc, and they represent the points at which | |
| 34 ## the array @var{vi} is interpolated. | |
| 35 ## | |
| 7001 | 36 ## If @var{x1}, @dots{}, @var{xn} are omitted, they are assumed to be |
| 6702 | 37 ## @code{x1 = 1 : size (@var{v}, 1)}, etc. If @var{m} is specified, then |
| 7001 | 38 ## the interpolation adds a point half way between each of the interpolation |
| 6702 | 39 ## points. This process is performed @var{m} times. If only @var{v} is |
| 40 ## specified, then @var{m} is assumed to be @code{1}. | |
| 41 ## | |
| 42 ## Method is one of: | |
| 43 ## | |
| 44 ## @table @asis | |
| 45 ## @item 'nearest' | |
| 46 ## Return the nearest neighbour. | |
| 47 ## @item 'linear' | |
| 48 ## Linear interpolation from nearest neighbours. | |
| 49 ## @item 'cubic' | |
| 50 ## Cubic interpolation from four nearest neighbours (not implemented yet). | |
| 51 ## @item 'spline' | |
| 52 ## Cubic spline interpolation--smooth first and second derivatives | |
| 53 ## throughout the curve. | |
| 54 ## @end table | |
| 55 ## | |
| 56 ## The default method is 'linear'. | |
| 57 ## | |
| 8116 | 58 ## If @var{extrapval} is the scalar value, use it to replace the values |
| 59 ## beyond the endpoints with that number. If @var{extrapval} is missing, | |
| 60 ## assume NA. | |
| 6702 | 61 ## @seealso{interp1, interp2, spline, ndgrid} |
| 62 ## @end deftypefn | |
| 63 | |
| 64 function vi = interpn (varargin) | |
| 65 | |
| 66 method = "linear"; | |
| 6742 | 67 extrapval = NA; |
| 6702 | 68 nargs = nargin; |
| 69 | |
| 70 if (nargin < 1) | |
| 71 print_usage (); | |
| 72 endif | |
| 73 | |
| 7208 | 74 if (ischar (varargin{end})) |
| 75 method = varargin{end}; | |
| 6702 | 76 nargs = nargs - 1; |
| 7208 | 77 elseif (ischar (varargin{end - 1})) |
| 78 if (! isnumeric (varargin{end}) || ! isscalar (varargin{end})) | |
| 6702 | 79 error ("extrapal is expected to be a numeric scalar"); |
| 80 endif | |
| 7208 | 81 method = varargin{end - 1}; |
| 8116 | 82 extrapval = varargin{end}; |
| 6702 | 83 nargs = nargs - 2; |
| 84 endif | |
| 85 | |
| 86 if (nargs < 3) | |
| 7208 | 87 v = varargin{1}; |
| 6702 | 88 m = 1; |
| 89 if (nargs == 2) | |
| 7208 | 90 m = varargin{2}; |
| 6702 | 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 | |
| 7208 | 103 elseif (! isvector (varargin{1}) && nargs == (ndims (varargin{1}) + 1)) |
| 104 v = varargin{1}; | |
| 6702 | 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 == | |
| 7208 | 113 (2 * ndims (varargin{ceil (nargs / 2)})) + 1) |
| 6702 | 114 nv = ceil (nargs / 2); |
| 7208 | 115 v = varargin{nv}; |
| 6702 | 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) = ":"; | |
| 6721 | 131 x{i} = x{i}(idx{:})(:); |
| 6702 | 132 endfor |
| 133 idx (1 : nd) = {1}; | |
| 134 idx (1) = ":"; | |
| 6721 | 135 x{1} = x{1}(idx{:})(:); |
| 6702 | 136 endif |
| 137 | |
| 7421 | 138 method = tolower (method); |
| 6702 | 139 |
| 140 if (strcmp (method, "linear")) | |
| 141 vi = __lin_interpn__ (x{:}, v, y{:}); | |
| 6742 | 142 vi (isna (vi)) = extrapval; |
| 6702 | 143 elseif (strcmp (method, "nearest")) |
| 144 yshape = size (y{1}); | |
| 145 yidx = cell (1, nd); | |
| 146 for i = 1 : nd | |
| 147 y{i} = y{i}(:); | |
|
7671
4fbaba9abec1
implement compiled binary lookup
Jaroslav Hajek <highegg@gmail.com>
parents:
7561
diff
changeset
|
148 yidx{i} = lookup (x{i}, y{i}, "lr"); |
| 6702 | 149 endfor |
| 150 idx = cell (1,nd); | |
| 151 for i = 1 : nd | |
| 7421 | 152 idx{i} = yidx{i} + (y{i} - x{i}(yidx{i}) > x{i}(yidx{i} + 1) - y{i}); |
| 6702 | 153 endfor |
| 154 vi = v (sub2ind (sz, idx{:})); | |
| 155 idx = zeros (prod(yshape),1); | |
| 156 for i = 1 : nd | |
| 157 idx |= y{i} < min (x{i}(:)) | y{i} > max (x{i}(:)); | |
| 158 endfor | |
| 159 vi(idx) = extrapval; | |
| 160 vi = reshape (vi, yshape); | |
| 6721 | 161 elseif (strcmp (method, "spline")) |
| 7423 | 162 if (any (! cellfun (@isvector, y))) |
| 163 for i = 2 : nd | |
| 164 if (! size_equal (y{1}, y{i})) | |
| 165 error ("dimensional mismatch"); | |
| 166 endif | |
| 167 idx (1 : nd) = {1}; | |
| 168 idx (i) = ":"; | |
| 169 y{i} = y{i}(idx{:}); | |
| 170 endfor | |
| 171 idx (1 : nd) = {1}; | |
| 172 idx (1) = ":"; | |
| 173 y{1} = y{1}(idx{:}); | |
| 174 endif | |
| 7421 | 175 |
| 6702 | 176 vi = __splinen__ (x, v, y, extrapval, "interpn"); |
| 7421 | 177 |
| 7423 | 178 if (size_equal (y{:})) |
| 179 ly = length (y{1}); | |
| 180 idx = cell (1, ly); | |
| 181 q = cell (1, nd); | |
| 182 for i = 1 : ly | |
| 183 q(:) = i; | |
| 184 idx {i} = q; | |
| 185 endfor | |
| 186 vi = vi (cellfun (@(x) sub2ind (size(vi), x{:}), idx)); | |
| 7424 | 187 vi = reshape (vi, size(y{1})); |
| 7423 | 188 endif |
| 6702 | 189 elseif (strcmp (method, "cubic")) |
| 190 error ("cubic interpolation not yet implemented"); | |
| 191 else | |
| 192 error ("unrecognized interpolation method"); | |
| 193 endif | |
| 194 | |
| 195 endfunction | |
| 196 | |
| 197 %!demo | |
| 198 %! A=[13,-1,12;5,4,3;1,6,2]; | |
| 199 %! x=[0,1,4]; y=[10,11,12]; | |
| 200 %! xi=linspace(min(x),max(x),17); | |
| 7421 | 201 %! AI=linspace(min(y),max(y),26)'; |
| 6702 | 202 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"linear").'); |
| 203 %! [x,y] = meshgrid(x,y); | |
| 204 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; | |
| 205 | |
| 206 %!demo | |
| 207 %! A=[13,-1,12;5,4,3;1,6,2]; | |
| 208 %! x=[0,1,4]; y=[10,11,12]; | |
| 209 %! xi=linspace(min(x),max(x),17); | |
| 210 %! yi=linspace(min(y),max(y),26)'; | |
| 211 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"nearest").'); | |
| 212 %! [x,y] = meshgrid(x,y); | |
| 213 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; | |
| 214 | |
| 215 %!#demo | |
| 216 %! A=[13,-1,12;5,4,3;1,6,2]; | |
| 217 %! x=[0,1,2]; y=[10,11,12]; | |
| 218 %! xi=linspace(min(x),max(x),17); | |
| 219 %! yi=linspace(min(y),max(y),26)'; | |
| 220 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"cubic").'); | |
| 221 %! [x,y] = meshgrid(x,y); | |
| 222 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; | |
| 223 | |
| 224 %!demo | |
| 225 %! A=[13,-1,12;5,4,3;1,6,2]; | |
| 226 %! x=[0,1,2]; y=[10,11,12]; | |
| 227 %! xi=linspace(min(x),max(x),17); | |
| 228 %! yi=linspace(min(y),max(y),26)'; | |
| 229 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"spline").'); | |
| 230 %! [x,y] = meshgrid(x,y); | |
| 231 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; | |
| 232 | |
| 6721 | 233 |
| 234 %!demo | |
| 235 %! x = y = z = -1:1; | |
| 236 %! f = @(x,y,z) x.^2 - y - z.^2; | |
| 237 %! [xx, yy, zz] = meshgrid (x, y, z); | |
| 238 %! v = f (xx,yy,zz); | |
| 239 %! xi = yi = zi = -1:0.1:1; | |
| 240 %! [xxi, yyi, zzi] = ndgrid (xi, yi, zi); | |
| 241 %! vi = interpn(x, y, z, v, xxi, yyi, zzi, 'spline'); | |
| 242 %! mesh (yi, zi, squeeze (vi(1,:,:))); | |
| 243 | |
| 7421 | 244 |
| 245 %!test | |
| 246 %! [x,y,z] = ndgrid(0:2); | |
| 247 %! f = x+y+z; | |
| 248 %! assert (interpn(x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5]), [1.5, 4.5]) | |
| 249 %! assert (interpn(x,y,z,f,[.51 1.51],[.51 1.51],[.51 1.51],'nearest'), [3, 6]) | |
| 250 %! assert (interpn(x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5],'spline'), [1.5, 4.5]) | |
| 251 %! assert (interpn(x,y,z,f,x,y,z), f) | |
| 252 %! assert (interpn(x,y,z,f,x,y,z,'nearest'), f) | |
| 253 %! assert (interpn(x,y,z,f,x,y,z,'spline'), f) | |
|
7561
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7424
diff
changeset
|
254 |
|
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7424
diff
changeset
|
255 %!test |
|
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7424
diff
changeset
|
256 %! [x,y,z] = ndgrid(0:2); |
|
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7424
diff
changeset
|
257 %! f = x.^2+y.^2+z.^2; |
|
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7424
diff
changeset
|
258 %! assert (interpn(x,y,-z,f,1.5,1.5,-1.5), 7.5) |