1% Create a sparse interpolation matrix.
2%
3% [H,out] = sparse_interp(mask,I)
4%
5% Create interpolation matrix from mask and fractional indexes I
6%
7% Input:
8%   mask: 0 invalid and 1 valid points (n-dimensional array)
9%   I: fractional indexes (2-dim array n by mi, where mi is the number of points to interpolate)
10% Ouput:
11%   H: sparse matrix with interpolation coefficients
12%   out: true if value outside of grid
13
14function [H,out] = sparse_interp(mask,I,iscyclic)
15
16if ndims(mask) ~= size(I,1) && ~(isvector(mask) && size(I,1) == 1)
17  error('sparse_interp: inconsistent arguments')
18end
19
20if isvector(mask)
21  sz = numel(mask);
22  mask = reshape(mask,[sz 1]);
23else
24  sz = size(mask);
25end
26
27% n dimension of the problem
28n = size(I,1);
29m = prod(sz);
30
31if nargin == 2
32  iscyclic = zeros(1,n);
33end
34
35
36% mi is the number of arbitrarly distributed observations
37mi = size(I,2);
38
39% handle cyclic dimensions
40for i = 1:n
41  if iscyclic(i)
42    % bring I(i,:) inside the interval [1 sz(i)+1[
43    % since the i-th dimension is cyclic
44
45    I(i,:) = mod(I(i,:)-1,sz(i))+1;
46  end
47end
48
49
50scale = [1 cumprod(sz(1:end-1))];
51
52% integer index
53ind = floor(I);
54inside = true(1,mi);
55
56for i = 1:n
57  if ~iscyclic(i)
58    % make a range check only for non-cyclic dimension
59
60    % handle border cases
61    p = find(I(i,:) == sz(i));
62    ind(i,p) = sz(i)-1;
63
64    inside = inside & 1 <= ind(i,:) & ind(i,:) < sz(i);
65  end
66end
67
68%whos ind
69
70% interpolation coefficient
71coeff(:,:,2) = I - ind;
72coeff(:,:,1) = 1 - coeff(:,:,2);
73
74% consider now only points inside
75iind = find(inside);
76coeff2 = coeff(:,iind,:);
77ind2 = ind(:,iind);
78mip = length(iind); % number of obs. inside
79
80si = repmat(1:mip,[2^n 1]);
81sj = ones(2^n, mip);
82ss = ones(2^n, mip);
83
84% loop over all corner of hypercube
85for i=1:2^n
86
87  % loop over all dimensions
88  for j=1:n
89    bit = bitget(i,j);
90
91    % the j-index of the i-th corner has the index ip
92    % this index ip is zero-based
93    ip = (ind2(j,:) + bit - 1);
94
95    % ip must be [0 and sz(j)-1] (zero-based)
96    % we know aleady that the point is inside the domain
97    % so, if it is outside this range then it is because of periodicity
98    ip = mod(ip,sz(j));
99
100    sj(i,:) = sj(i,:) + scale(j) * ip;
101    ss(i,:) = ss(i,:) .* coeff2(j,:,bit + 1);
102  end
103end
104
105% sj must refer to a valid point or its interpolation coefficient
106% must be zero
107
108inside(iind) = all(mask(sj) | ss == 0,1);
109H = sparse(iind(si(:)),sj(:),ss(:),mi,m);
110
111out = ~inside;
112% LocalWords:  interp
113
114% Copyright (C) 2004 Alexander Barth <a.barth@ulg.ac.be>
115%
116% This program is free software; you can redistribute it and/or modify it under
117% the terms of the GNU General Public License as published by the Free Software
118% Foundation; either version 2 of the License, or (at your option) any later
119% version.
120%
121% This program is distributed in the hope that it will be useful, but WITHOUT
122% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
123% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
124% details.
125%
126% You should have received a copy of the GNU General Public License along with
127% this program; if not, see <http://www.gnu.org/licenses/>.
128