1% Derive fractional indices on a separable grid. 2% 3% I = localize_separable_grid(xi,mask,x) 4% 5% Derive fractional indices where xi are the points to localize in the 6% separable grid x (every dimension in independent on other dimension). 7% The output I is an n-by-m array where n number of dimensions and m number of 8% observations 9 10function I = localize_separable_grid(xi,mask,x) 11 12x = cat_cell_array(x); 13xi = cat_cell_array(xi); 14 15% m is the number of arbitrarily distributed observations 16tmp = size(xi); 17mi = prod(tmp(1:end-1)); 18 19% sz is the size of the grid 20tmp = size(x); 21sz = tmp(1:end-1); 22 23if isvector(x) 24 n = 1; 25 mi = length(xi); 26 I = zeros(1,mi); 27 I(1,:) = interp1(x,1:length(x),xi); 28else 29 % n dimension of the problem 30 n = tmp(end); 31 m = prod(sz); 32 33 xi = reshape(xi,mi,n); 34 x = reshape(x,[m n]); 35 36 IJ = cell(n,1); 37 vi = cell(n,1); 38 X = cell(n,1); 39 XI = cell(n,1); 40 I = zeros(n,mi); 41 42 for i=1:n 43 vi{i} = 1:sz(i); 44 X{i} = reshape(x(:,i),sz); 45 XI{i} = xi(:,i); 46 end 47 48 [IJ{:}] = ndgrid(vi{:}); 49 50 51 for i=1:n 52 I(i,:) = interpn(X{:},IJ{i},XI{:}); 53 end 54 55end 56 57 58% handle rounding errors 59% snap to domain bounding box if difference does not exceeds tol 60tol = 50*eps; 61 62for i=1:n 63 % upper bound 64 ind = sz(i) < I(i,:) & I(i,:) <= sz(i) + tol; 65 I(i,ind) = sz(i); 66 67 % lower bound 68 ind = 1 < I(i,:) & I(i,:) <= 1 + tol; 69 I(i,ind) = 1; 70end 71 72 73 74 75% LocalWords: indices sz tol 76 77% Copyright (C) 2014 Alexander Barth <a.barth@ulg.ac.be> 78% 79% This program is free software; you can redistribute it and/or modify it under 80% the terms of the GNU General Public License as published by the Free Software 81% Foundation; either version 2 of the License, or (at your option) any later 82% version. 83% 84% This program is distributed in the hope that it will be useful, but WITHOUT 85% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or 86% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 87% details. 88% 89% You should have received a copy of the GNU General Public License along with 90% this program; if not, see <http://www.gnu.org/licenses/>. 91