Mercurial > hg > octave-lyh
annotate scripts/general/interpn.m @ 7671:4fbaba9abec1
implement compiled binary lookup
| author | Jaroslav Hajek <highegg@gmail.com> |
|---|---|
| date | Fri, 28 Mar 2008 15:53:09 -0400 |
| parents | a938cd7869b2 |
| children | 6c69f5cda32b |
| rev | line source |
|---|---|
| 6702 | 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 | |
| 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. | |
| 28 ## Each element of then @var{n}-dimensional array @var{v} represents a value | |
| 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 ## | |
| 58 ## If @var{extrap} is the string 'extrap', then extrapolate values beyond | |
| 59 ## the endpoints. If @var{extrap} is a number, replace values beyond the | |
| 6742 | 60 ## endpoints with that number. If @var{extrap} is missing, 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}; |
| 6702 | 82 nargs = nargs - 2; |
| 83 endif | |
| 84 | |
| 85 if (nargs < 3) | |
| 7208 | 86 v = varargin{1}; |
| 6702 | 87 m = 1; |
| 88 if (nargs == 2) | |
| 7208 | 89 m = varargin{2}; |
| 6702 | 90 if (! isnumeric (m) || ! isscalar (m) || floor (m) != m) |
| 91 error ("m is expected to be a integer scalar"); | |
| 92 endif | |
| 93 endif | |
| 94 sz = size (v); | |
| 95 nd = ndims (v); | |
| 96 x = cell (1, nd); | |
| 97 y = cell (1, nd); | |
| 98 for i = 1 : nd; | |
| 99 x{i} = 1 : sz(i); | |
| 100 y{i} = 1 : (1 / (2 ^ m)) : sz(i); | |
| 101 endfor | |
| 7208 | 102 elseif (! isvector (varargin{1}) && nargs == (ndims (varargin{1}) + 1)) |
| 103 v = varargin{1}; | |
| 6702 | 104 sz = size (v); |
| 105 nd = ndims (v); | |
| 106 x = cell (1, nd); | |
| 107 y = varargin (2 : nargs); | |
| 108 for i = 1 : nd; | |
| 109 x{i} = 1 : sz(i); | |
| 110 endfor | |
| 111 elseif (rem (nargs, 2) == 1 && nargs == | |
| 7208 | 112 (2 * ndims (varargin{ceil (nargs / 2)})) + 1) |
| 6702 | 113 nv = ceil (nargs / 2); |
| 7208 | 114 v = varargin{nv}; |
| 6702 | 115 sz = size (v); |
| 116 nd = ndims (v); | |
| 117 x = varargin (1 : (nv - 1)); | |
| 118 y = varargin ((nv + 1) : nargs); | |
| 119 else | |
| 120 error ("wrong number or incorrectly formatted input arguments"); | |
| 121 endif | |
| 122 | |
| 123 if (any (! cellfun (@isvector, x))) | |
| 124 for i = 2 : nd | |
| 125 if (! size_equal (x{1}, x{i}) || ! size_equal (x{i}, v)) | |
| 126 error ("dimensional mismatch"); | |
| 127 endif | |
| 128 idx (1 : nd) = {1}; | |
| 129 idx (i) = ":"; | |
| 6721 | 130 x{i} = x{i}(idx{:})(:); |
| 6702 | 131 endfor |
| 132 idx (1 : nd) = {1}; | |
| 133 idx (1) = ":"; | |
| 6721 | 134 x{1} = x{1}(idx{:})(:); |
| 6702 | 135 endif |
| 136 | |
| 7421 | 137 method = tolower (method); |
| 6702 | 138 |
| 139 if (strcmp (method, "linear")) | |
| 140 vi = __lin_interpn__ (x{:}, v, y{:}); | |
| 6742 | 141 vi (isna (vi)) = extrapval; |
| 6702 | 142 elseif (strcmp (method, "nearest")) |
| 143 yshape = size (y{1}); | |
| 144 yidx = cell (1, nd); | |
| 145 for i = 1 : nd | |
| 146 y{i} = y{i}(:); | |
|
7671
4fbaba9abec1
implement compiled binary lookup
Jaroslav Hajek <highegg@gmail.com>
parents:
7561
diff
changeset
|
147 yidx{i} = lookup (x{i}, y{i}, "lr"); |
| 6702 | 148 endfor |
| 149 idx = cell (1,nd); | |
| 150 for i = 1 : nd | |
| 7421 | 151 idx{i} = yidx{i} + (y{i} - x{i}(yidx{i}) > x{i}(yidx{i} + 1) - y{i}); |
| 6702 | 152 endfor |
| 153 vi = v (sub2ind (sz, idx{:})); | |
| 154 idx = zeros (prod(yshape),1); | |
| 155 for i = 1 : nd | |
| 156 idx |= y{i} < min (x{i}(:)) | y{i} > max (x{i}(:)); | |
| 157 endfor | |
| 158 vi(idx) = extrapval; | |
| 159 vi = reshape (vi, yshape); | |
| 6721 | 160 elseif (strcmp (method, "spline")) |
| 7423 | 161 if (any (! cellfun (@isvector, y))) |
| 162 for i = 2 : nd | |
| 163 if (! size_equal (y{1}, y{i})) | |
| 164 error ("dimensional mismatch"); | |
| 165 endif | |
| 166 idx (1 : nd) = {1}; | |
| 167 idx (i) = ":"; | |
| 168 y{i} = y{i}(idx{:}); | |
| 169 endfor | |
| 170 idx (1 : nd) = {1}; | |
| 171 idx (1) = ":"; | |
| 172 y{1} = y{1}(idx{:}); | |
| 173 endif | |
| 7421 | 174 |
| 6702 | 175 vi = __splinen__ (x, v, y, extrapval, "interpn"); |
| 7421 | 176 |
| 7423 | 177 if (size_equal (y{:})) |
| 178 ly = length (y{1}); | |
| 179 idx = cell (1, ly); | |
| 180 q = cell (1, nd); | |
| 181 for i = 1 : ly | |
| 182 q(:) = i; | |
| 183 idx {i} = q; | |
| 184 endfor | |
| 185 vi = vi (cellfun (@(x) sub2ind (size(vi), x{:}), idx)); | |
| 7424 | 186 vi = reshape (vi, size(y{1})); |
| 7423 | 187 endif |
| 6702 | 188 elseif (strcmp (method, "cubic")) |
| 189 error ("cubic interpolation not yet implemented"); | |
| 190 else | |
| 191 error ("unrecognized interpolation method"); | |
| 192 endif | |
| 193 | |
| 194 endfunction | |
| 195 | |
| 196 %!demo | |
| 197 %! A=[13,-1,12;5,4,3;1,6,2]; | |
| 198 %! x=[0,1,4]; y=[10,11,12]; | |
| 199 %! xi=linspace(min(x),max(x),17); | |
| 7421 | 200 %! AI=linspace(min(y),max(y),26)'; |
| 6702 | 201 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"linear").'); |
| 202 %! [x,y] = meshgrid(x,y); | |
| 203 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; | |
| 204 | |
| 205 %!demo | |
| 206 %! A=[13,-1,12;5,4,3;1,6,2]; | |
| 207 %! x=[0,1,4]; y=[10,11,12]; | |
| 208 %! xi=linspace(min(x),max(x),17); | |
| 209 %! yi=linspace(min(y),max(y),26)'; | |
| 210 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"nearest").'); | |
| 211 %! [x,y] = meshgrid(x,y); | |
| 212 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; | |
| 213 | |
| 214 %!#demo | |
| 215 %! A=[13,-1,12;5,4,3;1,6,2]; | |
| 216 %! x=[0,1,2]; y=[10,11,12]; | |
| 217 %! xi=linspace(min(x),max(x),17); | |
| 218 %! yi=linspace(min(y),max(y),26)'; | |
| 219 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"cubic").'); | |
| 220 %! [x,y] = meshgrid(x,y); | |
| 221 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; | |
| 222 | |
| 223 %!demo | |
| 224 %! A=[13,-1,12;5,4,3;1,6,2]; | |
| 225 %! x=[0,1,2]; y=[10,11,12]; | |
| 226 %! xi=linspace(min(x),max(x),17); | |
| 227 %! yi=linspace(min(y),max(y),26)'; | |
| 228 %! mesh(xi,yi,interpn(x,y,A.',xi,yi,"spline").'); | |
| 229 %! [x,y] = meshgrid(x,y); | |
| 230 %! hold on; plot3(x(:),y(:),A(:),"b*"); hold off; | |
| 231 | |
| 6721 | 232 |
| 233 %!demo | |
| 234 %! x = y = z = -1:1; | |
| 235 %! f = @(x,y,z) x.^2 - y - z.^2; | |
| 236 %! [xx, yy, zz] = meshgrid (x, y, z); | |
| 237 %! v = f (xx,yy,zz); | |
| 238 %! xi = yi = zi = -1:0.1:1; | |
| 239 %! [xxi, yyi, zzi] = ndgrid (xi, yi, zi); | |
| 240 %! vi = interpn(x, y, z, v, xxi, yyi, zzi, 'spline'); | |
| 241 %! mesh (yi, zi, squeeze (vi(1,:,:))); | |
| 242 | |
| 7421 | 243 |
| 244 %!test | |
| 245 %! [x,y,z] = ndgrid(0:2); | |
| 246 %! f = x+y+z; | |
| 247 %! assert (interpn(x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5]), [1.5, 4.5]) | |
| 248 %! assert (interpn(x,y,z,f,[.51 1.51],[.51 1.51],[.51 1.51],'nearest'), [3, 6]) | |
| 249 %! assert (interpn(x,y,z,f,[.5 1.5],[.5 1.5],[.5 1.5],'spline'), [1.5, 4.5]) | |
| 250 %! assert (interpn(x,y,z,f,x,y,z), f) | |
| 251 %! assert (interpn(x,y,z,f,x,y,z,'nearest'), f) | |
| 252 %! 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
|
253 |
|
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7424
diff
changeset
|
254 %!test |
|
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7424
diff
changeset
|
255 %! [x,y,z] = ndgrid(0:2); |
|
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7424
diff
changeset
|
256 %! f = x.^2+y.^2+z.^2; |
|
a938cd7869b2
__lin_interpn__.cc: handle decreasing coordinate values
Alexander Barth
parents:
7424
diff
changeset
|
257 %! assert (interpn(x,y,-z,f,1.5,1.5,-1.5), 7.5) |