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