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